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Abstract 


This  research  develops  a  feature-based  maximum  a  posteriori  (MAP)  classification 
system  and  applies  it  to  classify  several  common  pulse  compression  radar  and  communi¬ 
cation  modulations.  All  signal  parameters  are  treated  as  unknown  to  the  classifier  system 
except  SNR  and  the  signal  carrier  frequency.  The  features  are  derived  from  estimated  duty 
cycle,  cyclic  spectral  correlation,  and  cyclic  cumulants.  The  modulations  considered  in 
this  research  are  BPSK,  QPSK,  16-QAM,  64-QAM,  8-PSK,  and  16-PSK  communication 
modulations,  as  well  as  Barker5  coded,  Barker,  i  coded,  Barker511  coded,  Frank49  coded, 
PX49  coded,  and  LFM  pulse  compression  modulations.  Simulations  show  that  average  cor¬ 
rect  signal  modulation  type  classification  %C  >  90%  is  achieved  for  SNR  >  9dB,  average 
signal  modulation  family  classification  %C  >  90%  is  achieved  for  SNR  >  ldB,  and  an  av¬ 
erage  communication  versus  pulse  compression  radar  modulation  classification  %C  >  90% 
is  achieved  for  SNR  >  -4dB.  Also,  it  is  shown  that  the  classification  performance  using 
selected  input  features  is  sensitive  to  signal  bandwidth  but  not  to  carrier  frequency.  Mis¬ 
matched  bandwidth  between  training  and  testing  signals  caused  degraded  classification  of 
%C  ~  10%  -  14%  over  the  simulated  SNR  range. 


IV 


For  my  Family  and  Friends  who  listened  to  my  research  ramblings  and  helped  me  get 

through  the  stressful  days 


v 


Table  of  Contents 


Page 

Abstract .  iv 

Dedication .  v 

Table  of  Contents .  vi 

List  of  Figures . viii 

List  of  Tables .  x 

List  of  Acronyms .  xi 

I.  Introduction .  1 

1 . 1  Research  Motivation  and  Related  Research  .  1 

1 .2  Research  Goal .  2 

1.3  Research  Methodology .  3 

1 .4  Thesis  Organization .  4 

II.  Literature  Review  .  5 

2.1  Waveforms  Considered .  5 

2.1.1  Communication .  5 

2.1.2  Radar .  8 

2.2  Pattern  Recognition .  15 

2.2.1  Likelihood-Based  Tests .  15 

2.2.2  Feature  Based  Tests .  16 

2.3  Cyclostationarity .  16 

2.3.1  Theory .  18 

2.3.2  Cyclic  Autocorrelation  Function .  19 

2.3.3  Spectral  Correlation  Function . 22 

2.4  Estimating  the  Spectral  Correlation  Function . 23 

2.4.1  Temporal  Smoothing . 25 

2.4. 1.1  FFT  Accumulation  Method .  28 

2.4. 1.2  Strip  Spectral  Correlation  Algorithm . 30 

2.4.2  Frequency  Smoothing  .  31 

2.5  Cyclic  Cumulants .  31 


vi 


III.  Methodology 


Page 
.  36 


3.1  Simulating  Modulations .  36 

3.2  Simulating  SNR  with  AWGN . 41 

3.3  Extracting  Features . 43 

3.3.1  Duty  Cycle . 43 

3.3.2  Cyclic  Spectral  Correlation . 45 

3.3.3  Cyclic  Cumulants . 50 

3.4  Classifier  Training . 52 

3.5  Performance  Criteria . 55 

IV.  Results  and  Analysis . 58 

4.1  Simulation  Setup . 58 

4.2  Classifier  Performance  with  Ideal  Training  Data .  60 

4.2.1  Signal  Modulation  Type  Classification .  60 

4.2.2  Signal  Modulation  Family  Classification .  65 

4.2.3  Communication  vs.  Pulse  Compression  Radar  Modulation  Classi¬ 
fication  . 69 

4.3  Classifier  Bandwidth  Sensitivity . 73 

4.4  Classifier  Carrier  Frequency  Sensitivity  . 76 

V.  Conclusions . 77 

5.1  Summary . 77 

5.2  Impact . 78 

5.3  Recommendations  for  Future  Work . 79 

Bibliography  .  80 


vii 


List  of  Figures 


Figure 

2.1  Communication  Constellations . 

2.2  Pulse  Repetition  Interval . 

2.3  Nested  Barker4>5  Code  . 

2.4  Frequency  Spectrum  of  Frequency  Translates . 

2.5  SCF  Support  Region . 

2.6  Temporal  Smoothing . 

2.7  FAM  Estimate  Resolution . 

2.8  SSCA  Estimate  Resolution . 

2.9  Frequency  Smoothing . 

3.1  Waveform  Simulation  Process . 

3.2  MATLAB®  Generated  Pulse  Shaping  Filter  Properties  . 

3.3  MATLAB®  Generated  Pulse  Shaping  Filter  Applied  to  Simulated  BPSK  Signal 

3.4  Simulated  SNR  Scaling  Process  . 

3.5  Estimating  the  Duty  Cycle  in  Observation  Time  At  . 

3.6  Estimated  Duty  Cycles  Over  a  Range  of  SNR^g  with  95%  Confidence  Intervals 

3.7  Estimated  BPSK  SCF  at  SNR  =  20dB . 

3.8  Estimated  QPSK  SCF  at  SNR  =  20dB . 

3.9  Estimated  BPSK  SCF  at  SNR  =  -5 dB . 

3.10  Estimated  SCF  Feature  Ratio . 

3.11  Estimated  Cyclic  Cumulant  Spectrums  for  BPSK . 

3.12  Classifier  Training . 

3.13  Test  the  Classifier . 

3.14  Confusion  Matrix . 


Page 
.  7 
.  10 
.  14 
.  21 
.  23 
.  27 
.  29 
.  30 
.  32 
.  37 
.  39 
40 
.  41 
.  44 
45 
.  46 
.  47 
.  48 
.  50 
.  51 
.  54 
.  55 
.  56 


viii 


Figure  Page 

3.15  ROC  Curve  Examples . 57 

4. 1  Classifier  System’s  Average  Performance  for  12  Signal  Modulation  Types  with 

Ideal  Training  Data . 61 

4.2  Classifier  System’s  Modulation  Type  Classification  Performance  with  Ideal 

Training  Data .  62 

4.3  Classifier  System  ROCs  for  12  Modulation  Types  at  SNR  =  9dB . 63 

4.4  Classifier  System  ROCs  for  12  Modulation  Types  at  SNR  =  OdB . 64 

4.5  Classifier  System’s  Average  Performance  for  7  Modulation  Families  with  Ideal 

Training  Data .  66 

4.6  Classifier  System’s  Modulation  Family  Classification  Performance  with  Ideal 

Training  Data . 67 

4.7  Classifier  System  ROCs  for  7  Modulation  Families  at  SNR  =  OdB . 68 

4.8  Classifier  System’s  Average  Performance  for  Distinguishing  Communication 

from  Pulsed  Radar  Modulations  with  Ideal  Training  Data .  69 

4.9  Classifier  System’s  Pulsed  Radar  and  Communication  Modulation  Classifica¬ 
tion  Performance  with  Ideal  Training  Data . 70 

4.10  Classifier  System  ROCs  for  Communication  vs  Pulsed  Radar  Detection  at 

SNR  =  -5dB . 71 

4.1 1  Classifier  System’s  Performance  Sensitivity  to  Bandwidth . 73 

4.12  Classifier  System’s  Modulation  Type  Classification  Performance  with  Mis¬ 
matched  Bandwidth . 74 

4.13  Classifier  System  ROCs  for  12  Modulation  Types  at  SNR  =  8dB  with 

Mismatched  Bandwidth . 75 

4.14  Classifier  System’s  Performance  Sensitivity  to  Carrier  Frequency . 76 


IX 


List  of  Tables 


Table  Page 

2.1  Known  Barker  codes  [21] .  13 

2.2  Some  Frank  Code  Phase  Sequences  [21] .  14 

2.3  Some  P*  Code  Phase  Sequences  [21] .  14 

2.4  Cumulant  Partitions  for  n=4,  q= 2  33 

2.5  Cumulants .  34 

3.1  SCF  Classifier  Features . 49 

3.2  Cyclic  Cumulant  Features . 52 

3.3  Classifier  Features . 53 

4.1  Classifier  System’s  Confusion  Matrix  for  12  Modulation  Types  at  SNR  =  9dB  .  61 

4.2  Classifier  System’s  Confusion  Matrix  for  the  12  Modulation  Types  at  SNR  =  OdB  64 

4.3  Modulation  Families .  65 

4.4  Classifier  System’s  Confusion  Matrix  for  7  Modulation  Families  at  SNR=  OdB  68 

4.5  Radar  and  Communication  Waveforms . 70 

4.6  Classifier  System’s  Confusion  Matrix  for  Communication  vs  Pulsed  Radar 

Modulations  at  SNR  =  -ldB . 72 

4.7  Classifier  System’s  Confusion  Matrix  for  Communication  vs  Pulsed  Radar 

Modulations  at  SNR  =  -5dB . 72 

4.8  Classifier  System’s  Confusion  Matrix  for  12  Modulations  Types  at  SNR  =  8dB 

with  Mismatched  Bandwidths  . 75 


x 


List  of  Acronyms 


Acronym 

Definition 

ALRT 

average  likelihood  ratio  test 

ASK 

amplitude  shift  keying 

AWGN 

additive  white  gaussian  noise 

BPSK 

binary  phase  shift  keying 

CAF 

cyclic  autocorrelation  function 

CC 

cyclic  cumulant 

CTC 

cyclic  temporal  cumulant 

CTCF 

cyclic  temporal  cumulant  function 

DFT 

discrete  fourier  transform 

EW 

electronic  warfare 

FAM 

fast  fourier  transform  (FFT)  accumulation  method 

FFT 

fast  fourier  transform 

FSK 

frequency  shift  keying 

GLRT 

generalized  likelihood  ratio  test 

IF 

intermediate  frequency 

LFM 

linear  frequency  modulation 

MAP 

maximum  a  posteriori 

MATLAB® 

matrix  laboratory 

ML 

maximun  likelihood 

OFDM 

orthogonal  frequency  division  multiplexing 

PDF 

probability  density  function 

PMF 

probability  mass  function 

PRI 

pulse  repetition  interval 

xi 


Acronym 

Definition 

PSD 

power  spectral  density 

PSK 

phase  shift  keying 

QAM 

quadrature  amplitude  modulation 

QPSK 

quadrature  phase  shift  keying 

RADAR 

radio  detection  and  ranging 

RF 

radio  frequency 

ROC 

receiver  operating  characteristic 

SCF 

spectral  correlation  function 

SDR 

software  defined  radio 

SNR 

signal  to  noise  ratio 

SSCA 

strip  spectral  correlation  algorithm 

TCF 

temporal  cumulant  function 

TMF 

temporal  moment  function 

WSCS 

wide-sense  cyclo-stationary 

WSS 

wide- sense  stationary 

xii 


AUTOMATIC  MODULATION  CLASSIFICATION  OF  COMMON 


COMMUNICATION  AND  PULSE  COMPRESSION  RADAR 
WAVEFORMS  USING  CYCLIC  FEATURES 

I.  Introduction 

This  chapter  summarizes  the  research  presented  in  this  thesis.  Its  motivation  and  goals 
are  explained,  as  well  as  the  assumptions  used  to  limit  the  problem’s  scope.  Last, 
the  organization  of  information  and  results  presented  in  this  thesis  are  explained. 

1.1  Research  Motivation  and  Related  Research 

In  this  digital  age,  with  increasing  technology  and  decreasing  electronic  component 
size,  many  capabilities  are  being  integrated  into  single  complex  systems.  Also,  the 
ever  increasing  need  for  higher  data  rates  and  larger  bandwidths  in  the  electromagnetic 
spectrum  is  demanding  efficient,  adaptive  new  methods  to  utilize  the  licensed  and 
unlicensed  spectrums.  The  difficult  task  of  increasing  spectrum  usage  while  mitigating 
incurred  interference  between  independent  signals  can  benefit  from  automatic  modulation 
recognition  processes  applied  to  non-cooperative  signals  of  interest. 

Cognitive  radio  technology  with  software  defined  radios  (SDRs)  is  receiving  much 
research  interest  as  a  potential  solution  for  spectrum  management  problems  because  SDRs 
can  adaptively  change  critical  parameters  of  their  receive  and  transmit  operations  to  adjust 
to  current  channel  conditions.  Accurately  sensing  and  extracting  information  about  current 
spectrum  usage  is  a  key  process  for  a  cognitive  radio  system.  In  fact,  many  research  papers 
are  solely  focused  on  spectrum  sensing  techniques  for  cognitive  radios  [2,  18,  29].  The 
increasing  complexity  of  electromagnetic  environments  is  also  providing  new  challenges 
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for  electronic  warfare  (EW).  Spectrums  are  beginning  to  overlap  and  user  transmissions  are 
becoming  more  dynamic  in  time,  frequency,  and  modulation.  Improved  sensing  techniques 
of  the  electromagnetic  spectrum  is  key  for  future  communication  and  radar  systems  such 
as  cognitive  radios  and  cognitive  radars. 

Within  spectrum  sensing  research,  automatic  modulation  recognition  has  emerged 
as  an  important  process  in  cognitive  spectrum  management  and  EW  applications. 
Research  has  been  conducted  on  automatic  classification  of  both  digital  and  analog 
modulations  for  at  least  two  decades,  and  possible  applications  in  cognitive  radar 
and  communication  systems  include  threat  recognition  and  analysis,  communication 
interception/demodulation,  effective  adaptive  jammer  response,  and  communication/radar 
emitter  identification  [5,  23],  The  research  continues  to  trend  towards  larger  modulation 
sets  and  more  complicated  channel  environments  with  minimal  a  priori  signal  knowledge. 

In  [30],  the  feasibility  of  providing  automatic  modulation  recognition  as  an  integrative 
technology  for  radar  and  communication  signals  based  on  features  was  investigated, 
but  only  a  limited  set  of  modulation  types  were  simulated  and  varying  signal  to  noise 
ratio  (SNR)  analysis  was  not  provided.  [30]  is  the  only  research  found  that  addresses 
both  radar  and  communication  waveform  modulation  recognition.  This  area  remains 
relatively  unexplored  and  is  the  focus  of  this  research.  A  large  modulation  set  including 
both  pulse  compression  radar  and  communication  modulations  is  explored  for  modulation 
classification  with  minimal  knowledge  a  priori  of  critical  received  signal  parameters. 

1.2  Research  Goal 

The  goal  of  my  research  is  to  advance  the  application  of  modulation  classification 
presented  in  the  literature  by  developing  and  simulating  a  reliable  automatic  modulation 
recognition  system  capable  of  discerning  between  a  wide  range  of  non-cooperative  com¬ 
mon  pulse  compression  radar  and  communication  modulations.  Simulated  performance 
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and  limitations  of  the  developed  system  will  be  assessed  over  a  wide  range  of  received 
SNR  and  varying  received  signal  parameters. 

1.3  Research  Methodology 

First,  a  wide  set  of  communication  and  pulse  compression  radar  modulations  are 
simulated  with  varying  SNRs  by  adding  additive  white  gaussian  noise  (AWGN).  Then, 
promising  distinguishing  features  are  researched  and  chosen  for  use  in  a  classifier  system. 
The  research  is  directed  by  the  literature  which  documents  successful  feature-based 
classification  methods.  This  thesis  applies  these  research  findings  to  develop  and  simulate 
a  reliable  modulation  classification  system  for  both  common  communication  and  pulse 
compression  radar  modulations. 

In  [5],  a  survey  is  provided  of  prior  research  for  automatic  communication  modulation 
classification  techniques.  These  techniques  are  organized  by  statistical-based  and  feature- 
based  methods.  Although  statistical-based  techniques  are  theoretically  optimal,  they  are 
practically  inefficient  due  to  computational  complexity.  Feature-based  techniques  using 
cyclic  spectrum  features  and  cumulants  are  shown  to  have  performed  well  for  varying  sets 
of  communication  modulations  and  unknown  parameters.  These  same  parameters  were 
also  shown  to  perform  well  for  radar  waveform  modulation  recognition  in  [23],  and  [30] 
illustrated  that  the  estimated  duty  cycle  of  a  received  waveform  may  be  used  to  distinguish 
between  pulsed  radar  (linear  frequency  modulation  (LFM)  and  bi-phase  barker5  coded) 
and  conventional  communication  (AM,  FM,  ASK,  FSK,  BPSK,  QPSK)  signals  with  100% 
accuracy  for  SNR  greater  than  8dB. 

The  research  performed  in  this  thesis  is  focused  on  leveraging  signal  properties 
that  have  been  shown  to  be  successful  modulation  classification  features  to  develop  a 
versatile  classifier  system  capable  of  reliably  classifying  the  modulation  of  several  common 
communication  and  radar  modulated  waveforms.  These  signal  properties  include  signal 
duty  cycle,  cyclostationarity,  and  cyclic  cumulant  statistics  which  were  researched  for 
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classification  feature  selection  to  distinguish  between  binary  phase  shift  keying  (BPSK), 
quadrature  phase  shift  keying  (QPSK),  16-quadrature  amplitude  modulation  (QAM),  64- 
QAM,  8-phase  shift  keying  (PSK),  and  16-PSK  communication  signals  as  well  as  bi-phase 
Barker5  coded,  bi-phase  Barkern  coded,  bi-phase  Barker5>n  coded,  Frank49  coded,  Px49 
coded,  and  LFM  pulse  compression  radar  signals. 

1.4  Thesis  Organization 

Chapter  II  introduces  the  basic  theory  of  the  communication  and  radar  modulations 
considered,  and  describes  the  common  classification  methods  currently  utilized  for 
modulation  recognition.  It  then  summarizes  the  theory  found  in  literature  concerning 
cyclostationarity  and  various  algorithms  to  estimate  the  spectral  correlation  function  (SCF). 
Last,  the  topic  of  cyclic  cumulants  (CCs)  is  addressed. 

In  Chapter  III,  the  steps  taken  to  develop  the  modulation  classification  system  based 
on  the  theory  provided  in  Chapter  II  are  presented.  First,  the  process  of  simulating  the 
various  communication  and  radar  waveforms  is  explained  as  well  as  the  process  used  to 
simulate  the  received  SNR.  Then,  the  process  of  extracting  the  signal  features  researched 
in  Chapter  II  and  training  the  classifier  algorithm  is  explained.  Last,  the  criteria  used  to 
assess  the  developed  modulation  classification  system’s  performance  are  presented. 

Chapter  IV  provides  the  classifier’s  test  simulation  results  as  described  in  Chapter 
III.  Figures  for  probability  of  correct  classification  over  a  wide  SNR  range,  confusion 
matrices  for  SNRs  of  interest,  and  receiver  operating  characteristic  (ROC)  curves  for  SNRs 
of  interest  are  presented  for  multiple  test  simulations  with  varying  test  parameters.  These 
results  are  analyzed  and  compared  to  assess  the  classifier’s  performance. 

Finally,  Chapter  V  gives  a  summary  of  the  research  with  an  estimate  of  its  findings’ 
theoretical  and  operational  impact.  The  thesis  concludes  with  a  discussion  of  potential 
areas  for  continued  research  and  further  testing. 
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II.  Literature  Review 


This  chapter  provides  a  theoretical  background  of  concepts  used  in  this  research  as  well 
as  a  review  of  previous  work  published  in  the  area  of  interest.  Section  2. 1  covers 
the  various  communications  and  radar  waveforms  considered  in  the  research.  Section  2.2 
introduces  the  two  main  approaches  to  classification  and  pattern  recognition.  Section  2.3 
provides  the  development  of  cyclostationary  concepts  such  as  the  cyclic  autocorrelation 
function  (CAF)  and  spectral  correlation  function  (SCF).  These  concepts  are  extended  for 
practical  applications  by  introducing  various  methods  to  estimate  the  cyclic  spectrum  of 
signals  in  Section  2.4.  Last,  Section  2.5  provides  the  framework  for  higher-order  cyclic 
statistics  as  used  in  this  work. 

2.1  Waveforms  Considered 

This  research  includes  a  broad  range  of  common  communication  and  radar  waveforms 
for  modulation  recognition  analysis.  This  section  presents  the  fundamental  theory  for 
defining  each  modulation  type  and  provides  the  general  equations  that  represent  them. 
2.1.1  Communication. 

Digital  forms  of  communication  can  vary  envelope,  phase,  frequency,  or  any 
combination  of  these  to  relay  information  through  radio  frequency  (RF)  transmission. 
This  information  is  generally  encoded  and  represented  with  communication  symbols.  A 
modulation  scheme  utilizing  M  symbols  is  referred  to  as  M-ary.  The  simplest  modulation 
forms  only  modulate  in  one  domain  and  are  well  known  as  M-ary  amplitude  shift  keying 
(ASK),  M- ary  phase  shift  keying  (PSK),  and  M- ary  frequency  shift  keying  (FSK).  M-ary 
quadrature  amplitude  modulation  (QAM)  is  a  form  of  modulation  in  which  both  amplitude 
and  phase  are  varied  to  form  communication  symbols.  In  this  research,  binary  phase  shift 
keying  (BPSK),  quadrature  phase  shift  keying  (QPSK),  8-PSK,  16-PSK,  16-QAM,  and 
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64-QAM  communication  modulations  are  considered.  All  theory  in  this  section  is  derived 
from  information  found  in  [24,  26]. 

M- ary  ASK  transfers  information  through  its  amplitude  where  each  amplitude  level 
represents  a  communication  symbol.  Transmitted  ASK  symbols  at  a  carrier  frequency,  fc, 
fit  the  mathematical  form 


ASK  :  s(t )  =  Am  cos  (2 nfct)  0  <  t  <  Tsym  (2.1) 

where  Am  is  one  of  M  distinct  envelope  amplitudes  representing  communication  symbols. 
In  M-ary  PSK,  however,  amplitude  is  constant  because  information  is  transferred  through 
its  carrier  phase.  The  carrier  phase  may  have  one  of  M  values  in  any  symbol  period  Tsym 
given  by 

0m  =  ^  (m  -  1)  (2.2) 

where  m  =  1, 2,  •  •  •  ,  M.  Therefore,  the  modulated  M- PSK  waveform  at  a  carrier  frequency 
is 

M-PSK  :  sm(t )  =  A  cos  ^2 nfct  +  ^  _  l)j  (2.3) 

where  0  <  t  <  Tsym  and  m  =  1,2, •••  ,M.  From  Equation  (2.3)  we  can  calculate  the 
transmitted  symbols  for  BPSK,  QPSK,  8-PSK  and  16-PSK: 


BPSK  :  sm(t)  =A  cos  (2 nfct  +  inn) 


QPSK  :  sm(t )  =A  cos  ^2 nfct  +  ^  (m  -  1)  j 
8-PSK  :  sm(t)  =A  cos  (2nfct  +  ^  (m  -  l)j 


16-PSK  :  sm(t)  =A  cos  \2nfct  +  -(m  - 


m  =  1,2  (2.4a) 

m  =  1,2, 3,4  (2.4b) 

m  =  1, 2,  •  •  •  ,8  (2.4c) 

m=  1,2,  -,16  (2.4d) 


where  0  <  t  <  Tsym.  Utilizing  the  trig  identity  cos(A  +  B)  =  cos  A  cos  B  -  sin  A  sin  B, 


these  transmitted  signals  can  be  represented  in  quadrature  form  with  the  basis  functions 
(t)  =  cos(2n fct)  and  0i(O  =  sin(2n fct).  The  signal  constellations  are  shown  in  this  two 
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o 


o 


-A  A 

(a)  BPSK  Constellations 


(e)  16  QAM  Constellations 


(d)  16  PSK  Constellations 


Figure  2. 1 :  Communication  Constellations 


7 


dimensional  basis  function  space  in  Figure  2.1.  Also,  M-ary  PSK  signals  have  constant 
envelope  magnitudes  so  the  constellation  points  are  equally  spaced  on  a  circle  of  radius  A 
centered  at  the  origin.  The  constellations  for  BPSK,  QPSK,  8-PSK,  and  16-PSK  are  shown 
in  Figure  2.1a,  Figure  2.1b,  Figure  2.1c,  and  Figure  2.1d  respectively. 

Figure  2.1a  illustrates  that  BPSK  can  be  equal  to  2  level  antipodal  binary  ASK.  For 
the  case  that  its  two  phases  are  separated  by  180°,  Equation  (2.4a)  for  BPSK  is  equal 
to  Antipodal  2-ASK  from  Equation  (2.1).  For  example,  let  the  phase  take  the  values 
6  =  0  and  n  so  that  the  transmitted  BPSK  signal  is 

sm(t)  =  A  cos  (2nfct  +  mn )  m  =  1,2  (2.5a) 

=  ±A  cos  (2  nfct)  (2.5b) 

It  can  be  seen  that  although  the  phase  is  being  shifted,  the  amplitude  of  the  BPSK  signal 
envelope  can  take  the  two  values  ±A  as  in  ASK. 

M-ary  QAM  varies  both  its  carrier  phase  and  envelope  amplitude  to  represent  data 
symbols.  M-QAM  modulated  signals  can  be  defined  as 

M-QAM  :  sm(t)  =  Am<f>i(t )  +  Z?m02(O  ,  0  <  t  <  Tsym  ,  m  =  1, 2,  •  •  •  ,  M  (2.6) 

where  (p\(t)  =  cos (2 nfct),  (0  =  sin (2nfct),  Am  and  Bm  are  defined  as  Am  =  (2 am  -  1)  - 
^^M  and  Bm  =  (2b m  -  1)  -  Vm  with  am  and  bm  all  combinations  of  integers  in  the  set 
[l,  2,  •  •  •  ,  Vm],  For  16-QAM,  Am  and  Bm  may  have  values  [-3,  -1, 1, 3]  and  for  64-QAM, 
Am  and  Bm  may  have  values  [-7,  -5,  -3,  -1, 1, 3, 5, 7],  The  constellations  for  16-QAM  and 
64-QAM  are  therefore  square  lattices  instead  of  circular  and  are  shown  in  Figure  2.1e  and 
Figure  2.  If  respectively. 

2.1.2  Radar. 

RADAR  stands  for  radio  detection  and  ranging  (RADAR)  and  it  summarizes  the  two 
main  tasks  of  RADAR  systems.  That  is  to  detect  targets  and  determine  their  range  from 
the  RADAR  system  [21],  The  selection  of  a  radar  waveform  and  its  specifications  are 
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fundamental  to  the  performance  and  capabilities  of  a  radar  system.  Generally,  the  received 
signal  energy  determines  the  reliability  of  detection,  but  the  specifications  of  the  waveform 
are  responsible  for  the  accuracy,  resolution,  and  ambiguity  of  range  and  Doppler  (range 
rate)  of  the  target  [21].  Variables  that  may  be  manipulated  in  RADAR  waveforms  include: 
operating  frequency,  peak  power,  pulse  duration,  bandwidth,  pulse  repetition  interval  (PRI), 
modulation  type/coding,  and  polarization.  In  general,  continuous  wave  RADAR  has  very 
good  Doppler  sensitivities  but  weak  range  resolution.  Pulsed  RADAR  is  very  versatile 
and,  depending  on  design,  can  have  good  radar  resolution  in  Doppler  and  range  estimates 
to  provide  both  long  range  detection  and  adequate  resolution.  Pulsed  waveforms  have 
dominated  radar  design  [21].  Also,  due  to  desirable  correlation  properties,  these  waveforms 
are  very  similar  to  common  communication  modulations.  Therefore,  this  research  focuses 
on  recognizing  linear  frequency  modulation  (LFM),  Barker  coded,  Frank  coded,  and  Px 
coded  pulse  compression  radar  modulations. 

First,  general  radar  equations  for  the  simple,  pulsed  sinusoid  as  presented  in  [28]  are 
included  to  illustrate  the  important  increased  performance  realized  with  pulse  compression 
modulations.  Pulse  compression  modulations  utilize  many  modulation  schemes  common 
in  communications  signals.  A  fundamental  parameter  of  RF  transmission  is  wavelength. 
RF  wavelength  is  a  function  of  the  speed  of  light,  c  and  carrier  frequency  fc. 


where  the  speed  of  light  c  =  3  x  108  m/s.  In  the  most  simplistic  sense,  a  single  tone  RADAR 
pulse  is  [28] 

s  ( t )  =  Rect j  cos  (2nfct) ,  0  <  t  <  r  (2.7) 

and  the  range  between  a  monostatic  radar  system  and  a  target  is  given  by  the  range  equation 
[28] 

cTr 

R  =  ~f  (2.8) 
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where  Tr  is  the  round  trip  time  of  the  radar  pulse.  Alternatively,  the  maximum  unambiguous 
range  is  [28] 


Unambiguous  Range  :  Ru  = 


c  PRI 

2 


(2.9) 


where  the  PRI  calculated  as  the  time  between  RADAR  pulses.  Figure  2.2  shows  two  radar 
pulses  in  an  observation  interval  At.  The  range  resolution  and  accuracy  is  determined  by 
the  pulse’s  duration  r,  the  speed  of  light  c,  and  the  received  signal  to  noise  ratio  (SNR)  as 
in  Equation  (2.10)  and  Equation  (2.11)  respectively  [28], 
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(2.10) 


Range  Resolution  :  A R  « 


c  T 
~2 


Range  Accuracy  :  dR  ~  — - 

2  V2  SNR 


(2.11) 


Range  resolution  represents  the  distance  required  between  two  distinct  targets  for  the 
RADAR  system  to  reliably  distinguish  between  them.  Equation  (2.9),  Equation  (2.10), 
and  Equation  (2.11)  provide  information  for  characterizing  the  performance  of  a  RADAR 
system.  Performance  improvements  in  RADAR  systems  have  been  towards  greater  spatial 
resolution  capabilities  of  targets  with  noisy  backgrounds  [28]. 

The  duty  cycle  of  a  constant  amplitude  pulsed  signal  is  the  ratio  of  the  average  transmit 
power  over  the  PRI  and  the  peak  transmit  power  within  a  pulse  [28]. 


avg 


PRI 


(2.12) 


The  average  transmit  power  Pavg  is  the  instantaneous  transmitted  pulse  power’s  integral 
over  the  PRI  divided  by  the  PRI  and  the  peak  transmitted  power  Po  is  calculated  as  the 
transmitted  pulse  power’s  maximum  value  over  the  pulse  interval  r. 

-PRI 


avg 


J-f 

PRI  Jo 


p(t)  dt 


P0  =  max  (pit)  |5) 


(2.13a) 

(2.13b) 


with  the  instantaneous  power  of  the  transmitted  pulse  p  (i t )  =  fv  (t)\2. 

Range  rate,  or  Doppler,  is  how  the  RADAR  determines  target  velocity  relative  to  the 
RADAR  system.  The  RADAR  to  target  range  rate,  resolution,  and  accuracy  are  given  by 
[28] 

Range-Rate  :  Rdot  =  ^  (2.14a) 

Range-Rate  Resolution  :  A Rdot  =  (2.14b) 

2  r 

Range-Rate  Accuracy  :  5R.  = - ■=  (2.14c) 

°  2tVTSNR 
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There  is  a  trade-off  between  range  and  range-rate  resolution  and  accuracy  determined 
by  the  pulse  length  r.  A  long  pulse  width  is  desired  for  acute  Doppler  resolution  and 
accuracy  while  a  short  pulse  is  desired  for  fine  range  resolution  and  accuracy.  However, 
pulsed  radar  can  achieve  both  good  range  and  range-rate  resolutions  through  the  use  of 
pulse  compression  techniques.  The  pulse  compression  modulations  considered  in  this 
research  are  LFM  chirped,  Barker  coded,  Frank  coded,  and  Px  coded  waveforms. 

Pulse  compression  waveforms  allow  the  receiver  to  separate  targets  with  overlapping 
received  pulse  returns.  A  compression  filter  is  used  to  produce  a  narrow  or  compressed 
pulse  from  the  pulse  compression  modulated  received  signal.  The  duration  of  the  pulse  is 
therefore  reduced  in  the  receiver  and  results  in  a  better  range  resolution  than  was  expected 
from  the  transmitted  pulse  duration  [28],  Therefore,  pulse  compression  modulation  grants 
the  increased  Doppler  range  resolution  of  a  long-pulse  while  retaining  the  range  resolution 
of  a  narrow-pulse  through  received  echo  processing  [8]. 

LFM  was  the  first  and  still  is  a  widely  used  pulse  compression  method.  In  LFM,  the 
frequency  of  the  signal  is  swept  linearly  during  the  pulse’s  duration  r  over  a  bandwidth  W 
at  the  rate  j.  The  effective  time-bandwidth  product  of  LFM  is  W  x  r  and  contributes  to  the 
increased  range  resolution  of  a  LFM  pulse  over  a  simple  sinusoidal  pulse.  The  equation  for 
LFM  is  [21] 

LFM  :  s(t )  =  Rect  ^-jcos^r  t^fo  +  ,  0  <  t  <  r  (2.15) 

where  W  is  the  bandwidth  that  is  linearly  swept  during  the  pulse  duration  r  and  /o  is  the 
center  frequency.  Using  a  pulse  compression  receiver,  the  range  resolution  is  [21] 

A R  w  —  (2.16) 

2  W 

which  is  dependent  on  the  LFM’s  bandwidth  instead  of  its  pulse  duration  as  in 
Equation  (2.10). 

The  next  few  pulse  compression  methods  use  phase-coded  RADAR.  Instead  of 
linearly  sweeping  frequency  in  a  pulse  duration  r,  phase-coding  divides  the  pulse  into  M 
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Table  2.1:  Known  Barker  codes  [21] 


Code  Length 

Code 

2 

11  or  10 

3 

110 

4 

1110  or  1101 

5 

11101 

7 

1110010 

11 

11100010010 

13 

1111100110101 

sub-pulses  which  are  assigned  a  phase  value  according  to  a  specific  phase  code  sequence. 
To  maintain  consistent  notation  with  the  communication  waveforms,  the  sub-pulse  duration 
will  be  referred  to  as  Tsym  and  is  calculated  as  Tsym  =  [21], 

The  next  pulse  compression  method  uses  a  very  popular  and  common  family  of  codes 
known  as  Barker  codes.  Barker  codes  of  M  length  yield  a  max  peak-to-peak  sidelobe  ratio 
of  M.  There  are  only  nine  known  Barker  codes  [21],  all  listed  in  Table  2.1;  however, 
Barker  codes  can  be  nested  to  produce  larger,  sub-optimal  sequences  such  as  the  length  20 
Barker45  nested  code  as  shown  in  Figure  2.3.  Bi-phase  Barker  coded  RADAR  waveforms 
are  expressed  as  [21] 

Bi-phase  Barker  :  sm(t)  =  Rect  j  cos  (2 nfct  +  cmn)  ,  mTsym  <  t  <  (m  +  1)  Tsym  (2.17) 

where  cm  is  the  mth  value  of  a  known  Barker  code  listed  in  Table  2.1. 

Frank  and  PA  codes  apply  for  phase  sequences  of  perfect  square  length  M  =  L2  where 
sm  for  (1  <m<  M )  is  equal  to  s(i^\)L+i2  for  1  <  l\  <  L  and  1  <  /2  <  L  .  These  phase  codes 
produce  improved  range-rate  resolution  and  accuracy  over  Barker  phase  codes  [21],  Their 
sequences  are  calculated  from  [21] 

S(h-\)L+i2(t)  =  cos  (2nfct  +  (f>iuh)  (2.18) 
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Barker4  =  1101  Barker  5  =  11101 


1101 

1 

Y 

1101 

1 

y 

1101 

1 

Y 

1101 

0 

1 

NOT 

1101 

1 

Y 

Barker 4 5  = 

110  1 

110  1 

110  1 

0  0  10 

110  1 

Figure  2.3:  Example  of  nested  Barker4  5  code 


Table  2.2:  Some  Frank  Code  Phase  Sequences  [21] 


Code  Length 

Code 

1 

0 

4 

0,  0,  0,  7T 

9 

0000  —  —  0  —  — 

Table  2.3:  Some  Px  Code  Phase  Sequences  [21] 


Code  Length 

Code 

1 

0 

4 

7T  ~7T  ~7T  7T 

4’  4  ’  4  ’  4 

9 

1 

O 

O 

O 

where 


Frank  :  (f>tuh  =  2n(h  -  1  )(/2  -  1)  /L 


Px  •  <Ph,l2  — 


even 


(2.19a) 

(2.19b) 


¥  [§-<*]  [¥*-'■]•  Loia 
Frank  phase  codes  produce  linearly  stepped  linear  phase  segments  as  do  Px  codes 
except  Px  codes  have  their  zero  phase-rate  segment  terms  in  the  middle  of  the  pulse  instead 
of  at  the  beginning  [21].  Phase  values  for  the  first  three  square  Frank  and  Px  phase  codes 
calculated  from  Equation  (2.19)  are  given  in  Table  2.2  and  Table  2.3  respectively. 
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2.2  Pattern  Recognition 

Pattern  Recognition  has  become  a  very  useful  tool  with  applications  in  many 
areas  including  electronic  warfare  (EW)  and  Cognitive  software  defined  radio  (SDR). 
Pattern  recognition  research  for  selecting  and  extracting  features,  developing  classifier 
learning  algorithms,  and  evaluating  classifier  performance  is  still  prevalent  in  the  literature 
[1,  5,  7,  19].  For  most  applications,  there  are  two  main  methods  of  pattern  recognition 
that  are  being  used  for  modulation  classification:  likelihood-based  and  feature-based.  The 
likelihood-based  approaches  strive  to  minimize  false  classification  and  theoretically  can 
achieve  near  optimal  performance,  but  are  impractical  in  application  due  to  computational 
complexity.  Feature-based  methods  are  much  more  computationally  efficient  and  have  been 
shown  to  achieve  near  optimal  performance  in  the  Bayesian  sense  [5].  A  survey  of  current 
literature  addressing  both  methods  as  applied  to  communication  modulation  classification 
was  presented  in  [5]  and  an  example  of  feature-based  classification  for  radar  waveform 
classification  has  been  presented  in  [23]. 

2.2.1  Likelihood-Based  Tests. 

Likelihood  based  classification  methods  hinge  on  accurately  modeling  the  signal 
of  interest  and  all  other  ‘non-signal’  components  that  comprise  the  received  signal  ’s 
probability  distribution.  Decisions  are  made  by  comparing  likelihood  ratios  against  a 
threshold.  Among  likelihood-based  approaches,  two  ways  to  model  the  received  signal’s 
probability  distribution  are  the  average  likelihood  ratio  test  (ALRT)  and  generalized 
likelihood  ratio  test  (GLRT)  [5],  Depending  on  the  information  known  a  priori  about  the 
signals  being  discriminated,  either  the  ALRT  or  the  GLRT  is  used. 

The  ALRT  method  treats  received  unknown  variables  as  random  variables  with 
assumed  known  probability  density  functions  (PDFs),  but  the  GLRT  method  treats  the 
received  unknown  variables  as  deterministic  unknowns.  Therefore  the  GLRT  method  does 
not  make  any  assumptions  about  the  signal  or  the  channel  parameters.  The  final  decision  is 
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then  based  on  a  maximun  likelihood  (ML)  comparison  [5,  22],  For  a  binary  classification 
problem, 

£i[H]\r(t)]  H] 

/  777  7~  ^  Xj  ,  J  =  A  (ALRT),  G  (GLRT)  (2.20) 

L^ok  1 (*)\ 

where  Aj  is  a  threshold  and  the  method  used  to  compute  the  likelihood  functions  £  forms 
either  the  ALRT  or  GLRT  on  the  left  side.  H\  represents  decision  ‘1’,  H0  represents 
decision  ‘0’  in  this  binary  case,  and  r(t )  is  the  received  waveform 

2.2.2  Feature  Based  Tests. 

Feature-based  classification  methods  use  extracted  statistics,  or  features,  from  a 
received  signal  to  make  classification  decisions  based  on  the  reduced  data  set.  This 
reduced  data  set  is  called  a  feature  vector  and  is  represented  by  ip.  Some  examples 
of  discriminating  features  include  symbol  rates,  signal  magnitude  variance,  duty  cycle, 
instantaneous  frequency,  instantaneous  phase,  cumulants,  and  many  others.  Many  feature- 
based  methods  require  some  a  priori  knowledge  of  signal  parameters  in  order  to  accurately 
calculate  signal  features.  The  extracted  signal  features  are  then  used  for  decision  making. 
Decision  making  methods  are  usually  based  on  feature  PDFs,  or  feature  vector  distances 
from  calculated  class  feature  vector  means  [5]. 

In  literature,  cyclostationary-based  features  have  gained  popularity  as  potential 
features  for  modulation  recognition  because  they  are  insensitive  to  unknown  signal  and 
channel  parameters  and  preserve  signal  phase  information  [22].  In  [27],  the  received 
signal’s  fourth-order  two  conjugate  cumulants  were  used  as  features  to  discriminate 
between  BPSK,  4- ASK,  16-QAM  and  8-PSK  when  carrier  phases,  frequency  offsets,  and 
timing  offsets  were  unknown. 

2.3  Cyclostationarity 

A  stationary  random  process  is  one  where  all  its  joint  moments  are  non- varying  and 
all  its  functions’  expected  values  are  stationary,  wide-sense  stationary  (WSS)  is  a  weaker 
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form  of  stationarity,  which  requires  only  the  1st  and  2nd  order  statistics  to  be  stationary 
(not  vary  with  a  shift  in  the  time  origin).  Therefore,  a  WSS  random  process  has  a  mean 
( jj.x )  and  autocorrelation  (Rx)  that  satisfy  the  following  conditions  [11,  20]: 

E  [x(T)]  =/ix  ,  V  t 

Rx(t,  r)  =Rx(t)  =  E  x  (t  +  0  x*  (t  -  0  ,  V  t 

where  r  is  some  time  delay.  Both  statistics  are  independent  of  the  time  origin  ( t )  and  the 
auto-correlation  function  only  depends  on  the  time  difference  (r)  between  samples.  All 
stationary  random  processes  are  WSS,  but  not  all  WSS  processes  are  stationary  [20], 

Instead  of  non- varying  means  and  autocorrelations,  wide- sense  cyclo- stationary 
(WSCS)  random  processes  have  periodic  means  and  autocorrelations  [15].  Therefore,  for 
cyclostationary  random  processes,  the  mean  (jux )  and  autocorrelation  (Rx)  are  periodic  for 
some  period  T()  and  satisfy  the  following  conditions  [11]: 

E  [. x(t  +  T0)]  =Rx(t  +  T0)  =  /ux(t) ,  V  t 
Rx(t  +  T0,  r)  =Rx(t,  r)  =  E  [*  (f  +  0  x*  (f  -  0]  ,  V  t 
RF  waveforms  commonly  exhibit  cyclostationary  properties  due  to  common  operations 
such  as  modulating,  coding,  multiplexing,  and  sampling  which  induce  periodicities  in  the 
statistics  of  the  signals.  The  periodicities  in  autocorrelation  produce  spectral  correlations 
which  can  be  exploited  for  signal  processing  [10]. 

To  accurately  calculate  /ux(t)  and  Rx(t,r),  we  would  have  to  use  ensemble  averaging 
over  many  observations  of  a  single  process  and  have  knowledge  of  PDFs.  However,  if  time 
averaging  over  a  single  observation  is  equal  to  ensemble  averaging  over  many  observations, 
the  random  process  can  be  described  as  ergodic.  It  is  a  reasonable  assumption  for  most 
waveforms  used  in  communication  and  radar  applications  that  the  first  and  second-order 
statistics  within  the  transmitted  waveform  satisfy  the  ergodic  property  [26],  Therefore, 
to  avoid  an  unnecessary  probabilistic  discussion,  signals  in  this  paper  are  assumed  to  be 
ergodic  in  the  mean  and  autocorrelation  function.  This  allows  us  to  treat  the  temporal 
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average  as  equivalent  to  the  expected  value,  or  ensemble  average  [20],  In  this  thesis, 
temporal  averaging  with  respect  to  t  will  be  denoted  as  <•), . 


E  [•]  =  <•>, 


lim 

T  — >oo 


i  r\ 

T  J-T/2 


N 


(•)  dt 


lim - 

N^oo  2N  + 


i2> 


n=-N 


(2.21) 


2.3.1  Theory. 

In  order  to  derive  the  mathematical  representation  of  cyclostationarity  and,  in  turn, 
produce  the  SCF,  it  is  easiest  to  start  from  simple  frequency  analysis.  Time  limited  and 
periodic  signals  can  be  expanded  into  a  summation  of  weighted  sinusoids  known  as  a 
Fourier  Series 

+oo  +oo 

x(t)  =  Yj  Xnej%nt  =  Xnej2nnfot,  n  e  I  (2.22) 

n——oo  vi— — co 

where  /0,  the  inverse  of  the  period,  T(h  is  the  fundamental  frequency,  I  denotes  an 
integer  set,  and  the  coefficients  Xn  are  the  sinusoidal  component  weights  at  frequencies 
/  =  jr  =  nfo.  Therefore,  if  a  signal  has  a  non-zero  Fourier  Series,  it  has  the  additive 
sinusoidal  components  of  frequency  /  with  weights  Xn  given  by 

1  rT/2 

Xn  4  lim  -  x(t)e-j2n"  fydt.  (2.23) 

r->oo  T  J-t/2 

Let  us  assume  that  the  signal  x(t)  contains  a  finite  frequency  component  given  by 
a  cos  (27 rat),  where  a  is  the  frequency  magnitude  component  at  f  =  a.  Therefore,  the 
complex  Fourier  Series  coefficient  of  x(t)  at  frequency  a  may  be  represented  by  [10, 1 1, 15] 


1  rT/2 

Ma  =  lim  —  I  x(t)  e  j2m>Tdt 
t^>oo  T  J-t/2 


(x(t)e~j2nat)t. 


(2.24a) 

(2.24b) 


and  the  resulting  coefficient  Af®,  for  the  theoretical  x(t)  at  frequency  a  equals  \  a. 
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2.3.2  Cyclic  Autocorrelation  Function. 

Now  let’s  progress  to  a  signal  produced  by  the  lag-product  of  another  signal.  This 
quadratic  transformation  produces 

yit,T)  =  x{t+^jx*{t-^j  (2.25) 

where  (•)*  denotes  the  complex  conjugate  and  r  is  a  time  delay.  The  signal  y(t)  contains 
additive  sinusoidal  components  if  and  only  if 

My  (r)  =  (yit,  r)  e-**M)t  =  (x  (f  +  x*  (f  -  ,  (2.26) 

is  non-zero  for  any  frequency  a  4-  0. 

It  may  be  apparent  that  Equation  (2.26),  the  Fourier  Coefficients  of  the  lag-product 
(r),  is  a  generalized  formula  of  the  conventional  autocorrelation  function  of  x,  Rx(t).  It 
can  be  shown  that  in  the  special  case  where  a  =  0,  M\t=u  (r)  is  equivalent  to  the  conventional 
autocorrelation  function  Rx(t). 


M;=0(t)  =  (y(t,  r)  e-^0t)t  =  (y(t,  r))t  =  (x  (f  +  0  x*  (f  -  ^ 
RAr)  =  E  [x  (< +  T~) X-  (,  -  I)]  =  (x  (,  +  l) X-  (< -  I)) 


(2.27a) 

(2.27b) 


Therefore,  M“=0(r)  =  RJt)  =  (x  (7  +  ^jx*  (V  -  and  M"(t)  may  be  interpreted  as 
an  autocorrelation  function  of  x(t)  with  a  cyclic  weighting  factor  of  e~jlmn .  In  literature, 
M"(j)  is  commonly  expressed  as  the  cyclic  autocorrelation  function  (CAF)  and  is  written 
as  [10,  11,  13,  15] 


R"(t)  =  lim  ]- 

rT/2 

Rx{T)e~j2natdt 

(2.28a) 

T  — >oo  T 

J-T/2 

=  {RAt) 

e~j2nat\ 

(2.28b) 

=  (^x(t  + 

IH- 1)  ‘~*~l 

(2.28c) 

By  definition,  Equation  (2.28)  is  not  identically  zero  as  a  function  of  r  if  and  only  if  x(t) 
contains  second-order  periodicity  with  frequency  a  ^  0.  Therefore,  the  CAF  highlights  the 
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second-order  periodicities  with  frequency  a  in  the  signal  x(t).  Also,  Equation  (2.28)  has 
the  same  form  as  Equation  (2.24)  which  tells  us  that  Rx(t)  is  a  Fourier  coefficient  in  the 
Fourier  series  expansion  of  Rx(t)  [10]. 

oo 

Rx(t)  =  Yj  Rax(r)ej2nat  ,  a  =  j  ,  n  el  (2.29) 

n--oo 

Instead  of  an  autocorrelation  function  with  a  cyclic  weighting  factor,  the  CAF  can  also 
be  interpreted  as  a  conventional  cross-correlation  between  two  identical  signals  separated 
by  a  in  frequency.  Let  u(t)  and  v(t)  be  the  signal  x(t)  multiplied  by  e±j2n‘i'  which  shifts  the 


frequency  components  of  x(t)  by  +|  as  illustrated  in  Figure  2.4. 

u(t)  =x(t)e~j2n^  (2.30a) 

v(t)  =x(t)ej2n ^  (2.30b) 

The  Fourier  transforms  of  u(t)  and  v(t)  show  that  their  frequency  spectrums  are 

u(f)  =  r  [um  =  r  [me-^]  =  x  (/  + 1 )  (2.31a) 

V(f)  =  T[v(t )]  =r[x(t)e+^]  =  x(/-  |)  (2.31b) 


and  the  Wiener-Khinchin  relation  tells  us  that  the  Fourier  transforms  of  Ru(j)  and  Rv(t) 
give  us  their  Power  Spectral  Densities  (PSDs)  [3,  12,  20]. 

S„(/)=S,(/+!)  (2.32a) 

5v(/)=5,(/-|).  (2.32b) 

Defining  u(t)  and  v(t)  as  frequency  shifted  versions  of  x(t)  leads  us  to  an  important 
conceptual  understanding  of  the  CAF.  It  can  be  shown  that  the  conventional  cross- 
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Figure  2.4:  Frequency  spectrum  of  frequency  translates  u(t )  and  v(t)  of  x(t) 


correlation  of  u(t)  and  v(t)  equals  the  CAF  of  x(t). 


Ruvij)  =  E 

=  i^u 

Ruvir)  =Rax(T ) 


x[t+T^)e~JKaU+Tl2)  x  xy-l) 


jna(t-T/2 ) 


(2.33a) 

(2.33b) 

(2.33c) 

(2.33d) 


This  illustrates  the  interpretation  that  the  CAF  is  simply  a  temporal  cross-correlation 
between  frequency-shifted  versions  of  a  signal. 
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2.3.3  Spectral  Correlation  Function. 

According  to  the  Wiener-Khinchin  and  cyclic  Wiener-Khinchin  relations,  the  Fourier 
transform  of  the  autocorrelation  function  is  the  power  spectral  density  (PSD)  and  the 
Fourier  transform  of  the  CAF  is  the  SCF  [12,  13]. 


Sx(f) 
Sax(f ) 


Rx(T)e-WTdT  =  T[Rx(r)] 


Rax(T)e~j2nfTdT  =  T[Rax(T)\ 


(2.34a) 

(2.34b) 


The  SCF  is  represented  on  a  bi-frequency  plane  because  it  is  a  function  of  both  frequency, 
/,  and  cyclic  frequency,  a.  Just  as  the  conventional  autocorrelation  function  is  a  special 
case  of  the  CAF  for  when  a  =  0,  the  PSD  is  included  in  the  SCF  as  the  special  case  when 
a  =  0.  Remember  from  Equation  (2.33)  that  the  cross-correlation  of  u(t )  and  v(t)  equals 
the  CAF  of  x(t).  It  follows  that 


SaAf)  =  r  [i Rax(r )]  =  r  [. Ruv(t )]  =  Suv(f)  (2.35) 

where  Suv(f )  is  the  spectral  density  of  cross  correlation  between  u(t)  and  v(t)  at 
the  frequency  /  and  Sx(f)  is  the  spectral  density  of  correlation  between  the  spectral 
components  of  x(t)  at  /  -  |  and  /  +  |.  The  SCF  of  x(t)  is  the  Fourier  transform  of  the 
temporal  cross-correlation  between  frequency-shifted  versions  of  x(t). 

Suppose  that  x(t)  in  u(t)  and  v(t)  in  Equation  (2.30)  are  band-limited  with  a  double¬ 
sided  bandwidth  2 B.  The  SCF  region  of  support  for  a  band-limited  signal  is  illustrated  in 
Figure  2.5.  At  the  cyclic  frequency  of  a  =  0,  all  spectral  components  of  the  correlated 
frequency  translates  of  x(t)  overlap.  However,  for  the  cyclic  frequency  a  =  -B,  only 
spectral  components  from  -f  to  f  overlap  and  therefore  Sx=~B(f )  only  supports  the 
frequency  region  <  /  <  f .  The  frequency  translates  have  no  overlapping  spectral 
components  when  |cr|  >  26. 


In  [13]  and  [14],  the  SCF  for  analog  and  digital  modulated  signals  are  derived.  It 
is  shown  that  signals  with  the  same  power  spectral  densities  may  have  distinct  cyclic 
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spectrums.  Also,  the  cyclic  spectrum  is  shown  to  be  robust  to  additive  white  gaussian  noise 
(AWGN)  because  stationary  noise  has  no  cyclic  correlation.  Therefore,  distinguishing 
signal  features  may  be  extracted  from  the  SCF  and  can  be  used  for  robust  classification 
in  varying  noise  environments.  Techniques  for  estimating  the  SCF  from  sampled  data  are 
explored  in  Section  2.4. 

2.4  Estimating  the  Spectral  Correlation  Function 

The  theoretical  SCF  equations  presented  thus  far  deal  with  signals  of  infinite  time 
duration.  In  practice,  only  finite  time  observations  of  a  signal  are  available  for  analysis  and, 
as  such,  a  substantial  amount  of  work  has  been  done  to  modify  the  underlying  equations  to 
produce  efficient,  accurate  SCF  estimates.  In  general,  temporal  and  frequency  smoothing 
are  the  two  methods  used  to  produce  these  estimates.  Both  methods  derive  from  the  SCF 
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cyclic  periodogram  estimate  [9,  10,  25]. 


where 


(2.36) 

✓>oo 

Xj(t,  /)  =  I  aT(t  -  u)x  ( u )  e~'2nfudu 

(2.37a) 

w+r/2 

=  1  x(u)  e~'2lTtudu 

Jt-T/ 2 

(2.37b) 

is  the  finite  time  Fourier  transform  of  x(t )  with  aT  (t  -  u),  a  data  tapering  window  of  width 
T.  In  the  context  of  Spectral  Correlation,  XT(t,  f )  is  commonly  referred  to  in  literature  as  a 
complex  demodulate.  For  statistical  reliability,  and  a  reliable  estimate,  the  time-bandwidth 
product  should  be  much  greater  than  1  (At  x  A /  »  1)  [25].  The  cyclic  periodogram  in 
Equation  (2.36)  has  a  temporal  resolution  dictated  by  the  data  tapering  window  aT  ( t  -  u ) 
in  XT(t,f )  giving  At  =  T.  The  frequency  resolution  is  also  dictated  by  the  data  tapering 
window  size,  A /  «  j  ~  The  resulting  time-bandwidth  product  of  Equation  (2.36)  is 
AtAf  *  At±  *  1. 

Applying  time-smoothing  to  Equation  (2.36)  gives  the  time-smoothed  cyclic  peri¬ 
odogram 

co 

SaXT(t,f)At=  J  S"T(u,f)  •  h\t(t  -  u)  du  (2.38a) 

—  CO 
CO 

=  J'  XT[u,f+^jX^[u,f-^jhAt(t-u)du  (2.38b) 


where  the  new  time  resolution  is  defined  by  At,  the  width  of  the  sliding  data  tapering 
window  function  hAr(t  -  u).  To  maintain  statistical  reliability,  the  data  tapering  window 
function  should  have  a  width  At  »  ^  «  T  so  that  the  time  bandwidth  product  AtxAf  »  1. 
Applying  frequency-smoothing  to  Equation  (2.36)  gives  the  frequency-smoothed  cyclic 
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periodogram 


SaXT(t,f)  A/=  J  SaXT(t,v)-hA/(f-v)dv  (2.39a) 

— oo 
oo 

=  f  Xt({’v+ ^)h&f(f~v)dv  (2.39b) 

— oo 

where  the  new  frequency  resolution  is  defined  by  A/,  the  bandwidth  of  the  bandpass 
filter  hA/(f  -  v).  To  maintain  statistical  reliability,  the  bandwidth  of  hAj  (  f  -  v)  should  be 
A /  »  ^  ^  so  that  the  time  bandwidth  product  At  x  A /  »  1.  It  can  be  shown  that  both 

the  time-smoothed  cyclic  periodogram  and  the  frequency-smoothed  cyclic  periodogram 
approach  perfect  estimations  of  the  SCF  when  the  following  limits  are  applied  [10,  13] 

S“(/)=  lim  lim  S%{tJ,)At  (2.40a) 

T  — >co  At—>oo 

=  lim  lim  SaXT(t,f,)Af  (2.40b) 

Both  estimates  produce  a  cyclic  frequency  resolution  Aa  ~  2.  anci  it  follows  that  to 
maintain  reliable  estimates  A/  »  2.  »  Ao'.  Therefore,  the  SCF  estimate  must  have  finer 
resolution  in  cyclic  frequency  (a)  than  in  spectral  frequency  (/)  to  be  statistically  reliable. 

The  time  smoothing  and  frequency  smoothing  methods  are  generally  well  suited  for 
different  applications  of  SCF  estimation.  In  general,  variants  of  the  time-smoothed  cyclic 
periodogram  are  well  suited  for  efficient  estimation  over  the  entire  bi-frequency  plane, 
whereas,  variants  derived  from  the  frequency- smoothed  cyclic  periodogram  are  more  suited 
for  estimating  the  SCF  at  particular  cyclic  frequencies  [22,  25], 

2.4.1  Temporal  Smoothing. 

All  temporal  smoothing  algorithms  for  estimating  the  SCF  are  derived  from  the 
temporally  smoothed  cyclic  periodogram  given  in  Equation  (2.38).  Incorporating  the  data 
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where  ckq  =  ft  ~  fi  ,  fo  =  ,  r  is  a  dummy  variable,  and  XT  [r,  ft  ]  is  the  discrete  version 

of  Equation  (2.37)  given  by 

N'- 1 

XT  [r,  ft]  =  V  aT  [to]  x  [r  +  to]  e~j2nfk{r+m)Ts  (2.43a) 

m- 0 

N'-l 

=  YjaT  [m]  X  [r  +  to]  e-j2nk(r+m)/N'  (2.43b) 

m=0 

where  x[n]  =  x  (t)  |/=J!rs,  Ts  =  jg,  fk  -  jjf--,  and  T  =  N'T s.  The  temporal  resolution 
At  =  NTS  and  the  frequency  resolution  A /  =  which  produces  a  time  bandwidth 
product  ArA/  =  ^  and  cyclic  frequency  resolution  Aar  ~  ^  For  statistical  reliability 

»  N'  and  Aa  A /  [25]. 

Equation  (2.41),  Equation  (2.42),  and  Figure  2.6  show  that  the  time  smoothed  cyclic 
cross  periodogram  is  basically  a  correlation  between  the  spectral  components  of  x  [n]  over 
the  time  observation  of  At.  The  time  smoothing  is  done  by  allowing  a  data  tapering  window 
of  length  T  time  to  slide  over  the  total  signal  observation  At  time  or  equivalently  a  data 
tapering  window  of  length  N'  samples  to  slides  over  the  total  data  samples  of  length  N. 
Again,  the  window  of  size  N'  samples  or  T  time  should  be  much  smaller  than  the  total 
observation  length  of  N  samples  or  At  time  for  statistically  accurate  estimates  [22,  25]. 
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Figure  2.6:  Temporal  Smoothing  [25] 


Since  this  method  is  deemed  computationally  inefficient,  ways  to  improve  the 
computational  efficiency  of  the  time-smoothed  spectral  estimates  were  explored  in  [25]. 
One  method  to  improve  the  computational  efficiency  is  to  decimate  Equation  (2.43) 
by  L,  where  L  <  N',  giving  Xr  [rL,  fk] .  This  reduces  the  number  of  correlations  in 
Equation  (2.42)  by  a  factor  of  L  from  N  to  P  =  j.  Equivalently,  instead  of  calculating 
Equation  (2.43)  N  times,  then  decimating  to  P  values,  a  system  can  simply  calculate  the  P 
values  of  XT  [rL,  fk\  by  shifting  x  [ n ]  by  L  samples  each  computation.  A  decimation  factor 
of  L  =  ^  has  been  shown  to  be  a  good  choice  to  increase  computational  efficiency  and 
minimize  adverse  effects  from  cycle  leakage  and  cycle  aliasing  [4].  The  time  smoothing 
with  decimation  cyclic  periodogram  is  [25] 

n+P 

Sax°T  [nL,  f0]At  =  J]XT  [rL,  fk]  X*T  [rL,  f]  hAt  [n  -  r]  (2.44) 
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where 


N'-l 

Xj  [rL,  fk\  =  YjaT  N]  X  [rL  +  m ]  e-i2nk{rL+m)IN'  (2.45) 

m= 0 

Another  method  to  improve  the  cyclic  spectral  estimates  computational  efficiency  is 
to  multiply  both  sides  of  Equation  (2.42)  with  the  sinusoidal  factor  e~j2nqm/N.  This  shifts 
the  left  side  in  cyclic  frequency  by  ^  =  qAa  where  q  =  [0, 1,  •  •  •  ,  N  -  1]  and  fits  the  right 
side  into  the  form  of  an  /V-point  fast  fourier  transform  (FFT)  [25]. 

N 

SI,  [ n ,  fo\At  (e-^qm/N)  =  YJXr  [r, /,]  X*T  [r, f]  hAr  [n  -  r ]  (e~^qrlN)  (2.46a) 

r= 0 

SaxlT+qAa  foU  =  r  [XT  [r,  ft]  X*T  [r,  fk]  *  hAt  [r]]^  (2.46b) 

=  r  [Xt  [r,  ft]  X*T  [r,  fk]]N  r  [hAt  kl]w  (2.46c) 

where  ★  denotes  a  convolution,  the  notation  T  [-],v  denotes  an  A-point  DFT,  XT  [r,fk]  is 
computed  as  in  Equation  (2.43),  f0  =  ^  (§),  and  a0  =  fk  -  ft  = 

Utilizing  the  concepts  above,  [25]  presents  the  FFT  accumulation  method  (FAM)  and 
strip  spectral  correlation  algorithm  (SSCA)  as  computationally  efficient  time  smoothing 
algorithms  to  estimate  the  cyclic  spectrum. 

2.4.1. 1  FFT  Accumulation  Method. 

The  FAM  applies  both  decimation  and  FFTs  to  Equation  (2.42),  resulting  in 

p-i 

Sa£qAa  [nL,  f0]At  =J]XT  [rL,  /,]  X*T  [rL,  f]  hAt  [n  -  r]  e~j2*qr/p  (2.47a) 

r=0 

T-l 

=  [rL,  f]  X*T  [rL,  f]  hAt  [n  -  r]  e~^qrL,N  (2.47b) 

T=  0 

=T  [XT  [rL,  fk]  X*T  [rL,  f,]]p  T  [hAt  [r]]P  {2  Ale) 

where  XT  [rL,  fk]  is  defined  as  in  Equation  (2.45),  a0  =  o,  +  qAa,  L  is  the  decimation 
factor,  and  P  =  The  time  and  frequency  resolutions  are  Aa  =  =  At  =  A-  =  j,  and 

A  /  [q]  =  Aa  -  \q\  Aa.  Since  A  a  =  jf,  A  f  [q]  can  be  reduced  to  (l  -  \q\  jf.  Therefore, 
the  time-bandwidth  product  is  AtAf  =  At  (  Aa  -  \q\  Aa)  =  ^  -  \q\. 
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Figure  2.7  shows  the  support  region  for  the  FAM.  To  minimize  the  point  estimates  near 
the  top  and  bottom  of  the  channel-pair  regions,  where  q  is  large  and  the  time  bandwidth 
product  is  reduced  resulting  in  less  reliable  estimates,  only  the  estimates  within  the  region 
center  ±Aa/2  are  retained.  This  leaves  only  the  terms  corresponding  to 

A  a  ,  A  a 

- <  qAa  <  — 

2  2 


N  N 

IN’  1  2  N' 

PL  PL 


IN 


-<  q  < 


2N' 


-  1 
-  1 


Therefore,  there  are  missing  estimates  for  some  cyclic  frequencies,  a,  where  the  estimates 
are  less  reliable.  These  missing  estimates  may  contain  important  cyclic  features  and 
therefore,  the  FAM  is  not  advised  when  location  of  cyclic  features  is  unknown  a  priori. 
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Figure  2.8:  SSCA  Estimate  Resolution  [22,  25] 


2. 4.1. 2  Strip  Spectral  Correlation  Algorithm. 

The  second  temporally-smoothed  cyclic  spectral  estimation  algorithm  is  the  SSCA, 
which  allows  estimates  of  all  cyclic  frequencies.  In  this  algorithm,  the  complex 
demodulates  XT  [n,  fl\  directly  multiply  with  x*(n),  which  produces  estimates  along  the 
frequency-skewed  line  a  =  2fk-2f0.  This  algorithm  has  been  shown  to  give  highly  efficient 
estimates  of  the  SCF  over  the  entire  bi-frequency  plane,  but  sacrifices  fine  frequency 
resolution  [25].  The  SSCA  is  given  by 

sM,Aa  =J^XT  [r,  fk]  x *  [r]  hAt  [n  -  r ]  e~j2^r/N  (2.48a) 

L  JAr  r=0 

=T  [XT  [r, ./,]  x*  [r]]w  r  [hAt  [r]]w  (2.48b) 

where  no  =  fk  +  #An  and  /o  =  7  -  The  temporal  and  frequency  resolutions  are  At  -  j, 
A  /  =  4  =  A,  and  Aa  «  A  =  A  making  the  time-bandwidth  product  At  A/  =  A..  Like  the 
FAM  algorithm,  let  N  »  N'  to  produce  reliable  estimates. 
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2.4.2  Frequency  Smoothing. 

The  frequency- smoothed  cyclic  periodogram  equation  was  given  in  Equation  (2.39). 
In  [4]  and  [22],  Equation  (2.39)  was  extended  to  discrete  sampled  time-series. 


I 

2  1 


[n,fo]Af  =  Yj  Xt 


N^_ 

"  2 


n,fk  +  — 


n,  fi  +  — 


hAf  [ r ] 


(2.49) 


where  uo  -  fk~  fi  ■,  fo  -  ,  hAf  [r]  represents  the  response  of  some  bandpass  filter  with 

bandwidth  A/,  and  the  complex  demodulate 

N- 1 

xt  [n,  fk ]  =  J]aT  [ m ]  x[n  +  m ]  e-j2nfk{n+m)Ts  (2.50) 

m= 0 


is  now  calculated  from  N  samples  instead  of  N'  samples.  Figure  2.9  gives  a  graphical 
representation  of  Equation  (2.49).  The  temporal  and  frequency  resolutions  for  the 
frequency- smoothed  SCF  are  A  /  =  y  -  At  =  T  =  NTS,  and  A  a  «  The  time- 
bandwidth  product  is  then,  AtAf  =  N'  and  we  let  N'  »  1  for  statistical  reliability.  It 
is  apparent  from  Equation  (2.49)  that  there  is  a  trade-off  between  statistical  reliability 
and  spectral  resolution  [25].  To  achieve  highly  reliable  SCF  estimates,  a  large  amount 
of  frequency  smoothing  is  desired,  but  if  the  spectrum  has  narrow  spectral  features,  the 
amount  of  spectral  smoothing  should  be  minimized  [25], 


2.5  Cyclic  Cumulants 

Statistics  are  used  to  describe  and  characterize  the  behavior  of  processes.  Specifically, 
the  moments  and  cumulants  of  processes  are  very  useful  for  describing  behavior.  Since 
cumulant  functions  generally  can  not  be  computed  from  experimental  time-series  data, 
they  are  usually  estimated  from  knowledge  of  moment  functions,  which  can  be  computed 
from  experimental  data  [3],  Temporal  and  spectral  cumulants  are  shown  theoretically  to 
exhibit  the  property  of  signal  selectivity  in  [16].  This  is  the  ability  to  to  detect  or  estimate 
parameters  of  a  specific  signal  in  a  received  waveform  even  when  corrupted  by  noise  or 
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Average  in 
Frequency 


Figure  2.9:  Frequency  Smoothing 


interference.  This  property  was  verified  through  simulations  in  [27].  The  temporal  moment 
function  (TMF)  for  zero  time-lag  is  [3] 

Rx  if,  t  =  =  E  [x(tr(‘  (*W]  (2.51) 


and  is  used  to  compute  n-order,  ^-conjugate  moments.  It  is  apparent  that  the  autocorrelation 
defined  in  Equation  (2.27b)  is  a  specific  case  of  the  TMF  with  n  =  2  and  q  =  1  so 
Rx(t,T)2t i  =  E  [x(t)  x*(t)].  Using  the  moments,  cumulants  are  calculated  through  the 

moment  to  cumulant  formula,  also  known  as  the  temporal  cumulant  function  (TCF)  [16, 17] 

[  p  1 


Cx  it,  T\q  =  ^  (-ly  1  ip  -  1)!  Y\  Rx 

j= 1 


Ij 


(2.52) 


where  Pn  are  all  distinct  partitions  of  the  set  [1, 2,  •  •  •  , «],  p  is  the  number  of  elements  in 
each  partition,  and  Rx  ( t ,  r)n  q  is  the  n-order,  ^-conjugate  moment  corresponding  to  the  jth 
element  in  the  partition  [16].  It  has  been  shown  that  the  cyclic  cumulants  attain  maximum 
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Table  2.4:  n= 4,  q= 2  cumulant  partitions  where  (•)*  denotes  a  conjugate  and  ‘1’  and  ‘2’ 
were  generically  chosen  as  the  two  conjugated  terms. 


p 

n=  4,  q  =  2  Partitions 

Pn 

Cx  00,2  Partitions 

(-ir1  (p-D!n^o,o)n.,. 

j= 1 

1 

(102*,  3,4) 

Rx  0)4,2 

2 

(E,2*)  (3,4) 

-Rx  0)2,2  Rx  0)2,0 

2 

(1*,3)  (2*,  4) 

—Rx  0)2,1 

2 

(1*,4)  (2*,  3) 

—Rx  0)2,1 

2 

(1*,2*,  3)  (4) 

—Rx  0)3,2  Rx  0)i,o 

2 

(1*,2*,4)  (3) 

-Rx  0)3,2  Rx  0)i,o 

2 

(103,4)  (2*) 

—Rx  0)3,1  Rx  0)i,i 

2 

(203,4)  (1*) 

—Rx  0)3,1  Rx  0)i,i 

3 

(102*)  (3)  (4) 

2  Rx  0)2,2  Rx  (Oi,o 

3 

(103)  (2*)  (4) 

2  Rx  0)2,1  Rx  0)\,i  Rx  0)i,o 

3 

(104)  (2*)  (3) 

2  Rx  0)2,1  Rx  0)i,i  Rx  0)i,o 

3 

(2*,  3)  (1*)(4) 

2  Rx  0)2,1  Rx  0)i,i  Rx  0)i,o 

3 

(2*,  4)  (1*)(3) 

2  Rx  0)2,1  Rx  0)i,i  Rx  (Oi,o 

3 

(3,4)  (E)  (2*) 

2  Rx  0)2,0  Rx  0)\i 

4 

(E)  (2*)  (3)  (4) 

— 6Rx  0)i,i  Rx  0)\x) 

values  for  zero  delay  values,  r  =  0  [6].  Therefore,  all  n- order  moments  and  cumulants 
in  this  research  are  calculated  for  r  =  0  and  r  will  be  omitted  from  the  notation.  For 
example,  Cx  ( t ,  r)4>2  will  be  expressed  as  Cx  0)42  .  In  Table  2.4  an  example  for  calculating 
the  terms  for  Cx  (t)42  from  Equation  (2.52)  is  shown.  There  are  15  distinct  partitions  of  the 
set  [1*,  2*,  3, 4]  where  there  are  n  =  4  items  and  q  =  2  are  conjugated.  Item  ‘1’  and  item 
‘2’  in  the  set  were  generically  chosen  as  the  two  conjugated  terms,  but  any  combination 
of  two  may  be  chosen  as  long  as  the  selections  are  maintained  throughout  the  derivation. 
Summing  the  Cx  (t)42  partition  terms  in  Table  2.4  gives  Equation  (2.53). 

Cx  (^4,2  =  tfx  (04,2  —  \Rx  (02,0 1  —  2Rx  (02,i  —  2 Rx  0)3,2  Rx  (Oi,o  —  2/?x  0)3,1  Rx  (0l,l 

+  2 Rx  0)2,2  Rx  (Oi,o  +  %Rx  0)2,1  \Rx  (Oi,o|  +  Rx  0)2,0  Rx  (Oil  _  6  |/?x  (Oi,o|  (2.53) 
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Table  2.5:  Cumulants 


r 

n,q 

Equation 

6^2,0 

^2,0 

C2,i 

Ri,i 

C4,o 

R4,o 

C4,  i 

R41  -  3C2,oC2,i 

C4,2 

/?4,2-|C2,o|2-2C21 

Qi.o 

^6,0  _  1  5C2.oC4,o  -  15C2  0 

Cg,i 

^6,1  -  10C2>oC4>i  -  5C2j C4,o  -  15C2jiC|0 

Q,  2 

Rea  -  C;0C4,o  -  8C2>iC4>i  -  6C2,()C4i2  -  3q;0C220  -  12C2>oC2j 

Q,  3 

/?6,3  -  3C2  0^4,1  -  9C2>iC4>2  -  3C2>oC4  l  -  9C20C2,iC2;o  -  6C2 1 

C8,o 

Rs, 0  -  28C2>0C6>o  -  35C420  -  210C20  -  105C*0 

The  cumulant  equation  is  greatly  simplified  when  central  moments  are  used  instead 
of  raw  moments,  or  the  process  is  known  to  be  a  zero  mean  process,  fix  =  Rx  (t,  0)j  0  = 
Rx(t,  0)1;1  =  0.  In  practical  situations,  a  signal  can  be  made  a  zero  mean  process  by 
subtracting  the  mean  from  it.  Cx  (t,  0)4  2  reduces  to 

cx  {t\ 2  =  Rx  2  -  I  Rx  (02, o|2  -  2  Rx  (02,1  •  (2.54) 

A  list  of  the  zero-mean  cumulant  equations  derived  from  Equation  (2.52)  as  functions 
of  lower  order  moments  and  cumulants  are  shown  in  Table  2.5.  Owing  to  the  symmetrical 
signal  constellations  considered,  the  nth-order  moments  for  n  odd  are  zero  and  therefore, 
the  nth-order  cumulants  for  n  odd  are  also  zero  and  have  been  dropped  from  the  cumulant 
equations  in  Table  2.5  [5]. 

Much  like  the  CAF  is  found  by  Fourier  transforming  the  autocorrelation  function,  the 
cyclic  temporal  cumulant  function  (CTCF)  is  produced  by  Fourier  transforming  the  TCF 
[16] 

oo 

C$(t,  0\q  =  J  Cx(t,0\qe-^dt  (2.55) 

—co 

which  gives  the  TCF’s  frequency  components  at  frequency  y 3 .  The  nth-order,  ^-conjugate 
cycle  frequencies  (CFs)  of  interest  are  at  (5  =  (n-  2 q)  fc  [6].  Since  AWGN  is  a  stationary, 
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zero-mean  Gaussian  process,  its  cumulants  are  time  independent  and  non-zero  only  for  the 
second  order.  Therefore,  AWGN  does  not  have  any  contribution  to  the  higher-order  ( n  >  3) 
cyclic  cumulants  (CCs)  of  a  received  signal  r{t).  Last,  the  magnitude  of  the  nth-order, 
^-conjugate  CC  is  robust  to  the  carrier  phase  and  timing  offsets  [5,  6], 
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III.  Methodology 


This  chapter  outlines  the  work  that  led  to  the  development  of  a  modulation  recognition 
system.  The  modulation  recognition  system  is  feature-based  and  designed  to 
discriminate  between  BPSK,  QPSK,  16-QAM,  64-QAM,  8-PSK,  16-PSK,  Barker5, 
Barker 1 1 ,  Barker55,  Franl^g,  Px49,  and  LFM  modulations  using  features  derived  from 
theory  in  Chapter  II.  All  simulations  were  done  in  discrete-time  with  matrix  laboratory 
(MATLAB®),  therefore  all  equations  will  be  presented  for  discrete-time. 

Section  3.1  describes  the  process  used  to  simulate  the  waveforms  and  Section  3.2 
describes  the  process  of  introducing  AWGN  to  the  waveforms  to  simulate  received  SNR. 
Section  3.3  highlights  how  the  features  were  estimated,  Section  3.4  explains  the  classifier 
supervised  training  process,  and  Section  3.5  gives  the  metrics  used  to  assess  the  classifier 
performance. 

3.1  Simulating  Modulations 

This  section  describes  the  process  used  to  simulate  the  waveforms  being  considered 
in  this  research.  The  process  used  to  simulate  the  waveforms  in  MATLAB®  is  shown  in 
Figure  3.1. 

Equation  (3.1)  is  the  discrete  version  of  Equation  (2.4)  and  is  used  to  simulate  the 
discrete  symbols  for  BPSK,  QPSK,  8-PSK,  and  16-PSK 


BPSK  :  sm  [n]  =A  cos  (2 nfcnTs  +  n(m-  1)) 

m  -  1,2 

(3.1a) 

QPSK  :  sm  [n\  =A  cos  ^2nfcnTs  +  ^  (m  -  l)j 

m  =  1,2, 3, 4 

(3.1b) 

8-PSK  :  sm  [ n]  =A  cos  ^2nfcnTs  +  ^  (m  -  1)  j 

m  =  1, 2,  •  •  •  ,8 

(3.1c) 

16-PSK  :  sm  [n\  =A  cos  [lnfcnTs  +  ^  (m  -  1)  j 

SO 

<N 

II 

s 

(3. Id) 
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e(-j2nfcnTs ) 


Simulated  Waveform 

s[n] 


Figure  3.1:  Waveform  Simulation  Process 


where  each  symbol  is  defined  on  the  discrete  interval  0  <  nTs  <  Tsym,  Tsym  is  the  symbol 
period,  and  Ts  is  the  discrete  sampling  period  related  to  sampling  frequency  by  =  7-. 
Equation  (3.2)  is  the  discrete  form  of  Equation  (2.6)  and  is  used  to  generate  the  16-QAM 
and  64-QAM  symbols 


M-QAM  :  sm  [ n ]  =  Am(/)i(nTs)  +  Bm(p2(nTs )  (3.2) 

0  <  nTs  <  Tsym  m=  1, 2,  •  •  •  ,  M 

where  <p\(t)  =  cos (2 nfcnTs),  (p2(t)  =  sin (2 nfcnTs),  Am  and  Bm  are  defined  as  Am  = 
(2 am  -  1)  -  VI?  and  Bm  =  (2b m  -  1)  -  ^^M  with  am  and  bm  all  combinations  of  integers 
in  the  set  Jl,  2,  •  •  •  ,  VMj.  For  16-QAM,  Am  and  Bm  may  have  values  [-3,-1, 1,3]  and  for 
64-QAM,  Am  and  Bm  may  have  values  [-7,  -5,  -3,  -1, 1, 3, 5, 7]. 

Since  the  information  symbols  in  communication  waveforms  can  be  modeled  as 
random,  each  generated  communication  symbol  is  uniformly  randomly  selected.  This 
process  was  simulated  by  using  the  MATLAB®  command  ‘randi’  which  uniformly 
randomly  selects  a  value  from  a  given  set.  The  resulting  communication  waveform,  .v  [n], 
consists  of  a  stream  of  random  symbols  sm  [n].  This  produces  pseudo  random  streams  of 
each  modulation  where  all  symbols  are  equally  likely  to  occur  every  symbol  period.  Also, 
the  signal’s  modulated  and  transmitted  information  bandwidth  is  W  =  2B  =  7==-. 

-L  sym 
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Equation  (3.3)  is  used  for  simulating  an  up-chirp  pulse  with  bandwidth,  B  centered  at 
carrier  frequency  fc  with  duration  r. 

LFM  :  s  [n\  =  Rect  ^  j  cos  ^ 2n  nTs  |/c  +  ^-nTs  jj  ,  0  <  nTs  <  r  (3.3) 

Equation  (3.4)  is  used  to  simulate  Barker5,  Barker 1 1 ,  and  Barker511  codes.  Barker5  and 
Barkern  use  the  codes  given  in  Table  2.1  corresponding  to  length  5  and  11  respectively, 
but  Barkei'5  i  i  uses  the  code  generated  by  ‘nesting’  the  length  5  code  within  a  length  11 
code  sequence  similar  to  the  Barker4  5  example  in  Figure  2.3. 

(nTs  \ 

Barker  :  sm  \ri\  =  Recty - jcos  (2 nfcnTs  +  cmn )  ,  mTsym  <  nTs  <  (m  +  1)  Tsym  (3.4) 

Equation  (3.5)  is  the  discrete  equation  for  Frank  and  Px  coded  radar  pulses  with  the  phases 
defined  as  in  Equation  (2.19). 


S(h-i)L+h  [n\  =  cos  (2nfnTs  +  (pluh) 


(3.5) 


The  radar  pulse  compression  modulations  are  considered  deterministic  not  random.  Their 
value  during  a  symbol  period  is  predetermined  by  the  specific  code  sequence  corresponding 
to  the  pulse  compression  format/type. 


After  the  waveform  symbols  are  generated,  they  are  filtered  using  a  raised  cosine  pulse 
shaping  filter  in  MATLAB®  with  50%  excess  bandwidth  and  a  roll-off  factor  /3  =  0.4.  This 
filter  was  used  to  simulate  the  pulse  shaping  filter  in  a  transmitter  and  to  band-limit  the 
transmitted  simulated  waveform.  The  impulse  response  of  a  raised  cosine  filter  is  given  by 

,(-AnLA 

\  1  sym  ) 


h  [ n\  =  sine 


nT , 


sym 


COS 


1  -  (2i^)2 

\  4  sym  ) 


(3.6) 


This  raised  cosine  filter  was  generated  in  MATLAB®  using  the  ‘firreos’  command  with  the 
options:  filter  order  =  10,  cutoff  frequency  =  1.5  x  B  where  B  =  fSym  =  jr—  ,  and  roll¬ 
off  factor  p  =  0.4.  It  is  applied  to  the  simulated  modulations  using  the  MATLAB®  ‘filtfilt’ 
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(a)  Impulse  Response  of  Pulse  Shaping  Filter  Used 


(b)  Frequency  Response  of  Pulse  Shaping  Filter 


Figure  3.2:  Pulse  shaping  filter  properties  using  MATLAB®  ‘firrcos’  command  with  order 
=  10,  cutoff  frequency  =  0.15/s,  and  roll-off  factor  [3  =  0.4  applied  with  the  ‘filtfilt’ 
command . 

command  which  effectively  squares  the  filter  response  by  applying  the  filter  twice  to  negate 
phase  distortion.  Thus,  the  pulse  shaping  filter  used  in  this  research  is  | h  [n]|2.  The  impulse 
and  frequency  response  for  this  MATLAB®  generated  pulse  shaping  filter  with  B  =  0.1  fs. 
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(a)  Temporal  representation  of  BPSK  filtered  by  pulse  shaping  filter.  (Red 
-  Filtered,  Blue  -  Unfiltered) 


(b)  Frequency  representation  of  BPSK  filtered  by  pulse  shaping  filter 
(Red  -  Filtered,  Blue  -  Unfiltered) 

Figure  3.3:  MATLAB®  generated  pulse  shaping  filter  from  Figure  3.2  applied  to  simulated 
BPSK  Signal  with  B  =  jr—  =  0.1/*. 


applied  with  the  ‘filtfilt’  command,  is  shown  in  Figure  3.2.  In  this  research,  bandwidth 
is  referred  to  as  both  B  and  W,  where  B  is  the  baseband  bandwidth  and  W  =  2B  is  the 
transmission  bandwidth. 
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N Noise  IV-] 


Figure  3.4:  Simulated  SNR  Scaling  Process 


Furthermore,  the  temporal  and  frequency  responses  from  passing  a  BPSK  signal  with 
W  =  2B  =  2 fs  through  the  pulse  shaping  filter  shown  in  Figure  3.2  is  shown  in  Figure  3.3. 
After  pulse  shaping,  the  signal  is  then  mixed  up  to  a  carrier  frequency  fc  by  multiplying  the 
baseband  signal  by  e{j2n-fc  nTs\ 

3.2  Simulating  SNR  with  AWGN 

After  a  waveform  is  simulated,  filtered,  and  upconverted  to  a  carrier  frequency  for 
transmission,  channel  effects  are  simulated  by  introducing  randomly  generated  AWGN  to 
simulate  a  specific  SNR.  The  SNR  is  simulated  through  the  process  shown  in  Figure  3.4. 

First  complex  AWGN,  NNoise  [n],  is  simulated.  AWGN  has  a  normal  Gaussian 
distribution.  This  was  realized  in  MATLAB®  by  generating  streams  of  random  real  and 
complex  values  using  the  ‘randn’  command.  The  simulated  waveform’s  power  and  the 
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simulated  AWGN’s  power  is  computed  using  Equation  (3.7) 


f W  =  l  £  K'*'2  =  <IOI2>„  =  [lOI2] 

m= 0 


and  the  ratio  of  signal  power,  P signal-,  and  the  noise  power,  P Noise,  is  the  SNR. 


SNR  = 


Psi 


S  ignal 


s[n] 


(3.7) 


(3.8) 


P  Noise  F,NNoi:se\n\ 

To  achieve  a  desired  SNR,  the  simulated  waveform,  s  [n\,  is  scaled  by  a  constant  scale 
factor,  scale,  found  through  the  relationship 


SNR 


Desired  ~ 


P S  ignal Desire(i 

P 

1  Noise 
P  s[n]xscale 

^ N Noise  \n\ 

[n]  x  scale^j 

(l NNoise  [n] \2)n 

x  (scale)2 


N Noise  In] 


(3.9a) 

(3.9b) 

(3.9c) 

(3.9d) 


Rearranging  and  solving  for  scale  in  the  above  equation  gives 

scale  =  -x  I  SNR 


Noise 


Desired " 


(3.10) 


P 

*  Signal 

so  that  multiplying  5  [n]  by  the  scale  factor  gives  us  the  desired  SNR. 

Therefore,  the  simulated  received  signal  r[n]  with  a  specific  SNR  is  produced  by 
multiplying  the  signal  (s  [n])  by  the  calculated  scale  factor  from  Equation  (3.10)  and  adding 
this  product  to  the  simulated  noise,  NNoise  [«]. 


r[n\  =  s  [n\  x  scale  +  NNoise  [n] 


(3.11) 


The  desired  SNRs  in  this  work  are  expressed  in  a  logarithmic  decibel  scale.  To  convert 
between  linear  and  decibel  scales  utilize  Equation  (3.12) 


SNRjB  =  101og10(SNR) 


SNR(/B 

SNR  =  10~ 


(3.12a) 

(3.12b) 
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3.3  Extracting  Features 

Once  the  received  signal  was  simulated,  feature  analysis  was  done  to  extract  useful 
features  for  use  in  classification.  The  following  subsections  describe  the  features  extracted 
from  the  simulated,  received  waveforms. 

3.3.1  Duty  Cycle. 

An  estimated  version  of  duty  cycle,  6C  is  used  to  determine  whether  the  received  signal 
is  present  during  the  whole  observation  period  or  is  a  pulse.  To  do  this,  |r  [n]|2  is  smoothed 
with  a  moving  average  filter  to  produce 


n+Ns  f~  1 


\r[n\ | 


smoothed 


Wr  Z  ' ^ 


(3.13) 


where  Ns  f  is  the  smoothing  factor  and  represents  the  number  of  samples  averaged  for  each 
smoothed  value.  This  research  uses  NSf  =  4  x  arbitrarily  chosen  based  on  simulation 
results.  6C  is  estimated  with  Equation  (2.12)  where  an  estimated  Pq  is  used  instead  of 
Equation  (2.13b)  and  Pavg  is  the  average  power  in  the  observed  time  interval  At.  P0  was 
chosen  to  be  estimated  as 


A)  =  max,,  [p  [n]smoothed]  =  maxn  [\r  [n] \]moothed\  ,  0  <  nts  <  At  (3.14) 

When  the  transmit  pulse  duration  r  is  not  known,  8C  provides  estimated  duty  cycles  that 
are  statistically  different  for  pulsed  and  continuous  waveforms  as  shown  in  Figure  3.5  . 
Therefore,  8C  is  calculated  as  ratio  between  the  average  power  in  an  observed  time  interval 
and  the  estimated  P0  from  Equation  (3.14) 

K  =  (3.15) 

P« 


The  received  waveform’s  estimated  duty  cycles  are  shown  in  Figure  3.6  for  various 
SNRs.  The  range  bars  show  the  95%  confidence  interval  for  the  estimated  duty  cycle  at  each 
SNRjg  based  on  3,000  simulated  observations  per  waveform  per  SNR.  It  can  be  seen  that 
the  estimated  duty  cycle  is  affected  by  noise  and  all  the  waveforms  tend  to  have  statistically 
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Figure  3.5:  Estimating  the  duty  cycle  of  an  arbitrary  pulse  with  SNR  =  20  dB  over  an 

observation  time  At  by  using  Po  instead  of  Pq  when  r  is  unknown.  6C  =  =  0.14  and 

§  _  ao7  =  o  13 
°c  0.53  u-i~5 


identical  estimated  duty  cycles  below  -2  SNRrfB.  It  can  also  be  seen  that  the  simulated 
communication  waveforms  have  statistically  greater  duty  cycles  than  the  simulated  pulse 
compression  radar  waveforms  at  SNRrfg  >  2.  These  results  are  similar  to  those  presented 
in  [30],  which  used  a  duty  cycle  threshold  of  0.4  to  separate  pulsed  radar  waveforms  from 
conventional  communication  waveforms.  Therefore,  6C  is  chosen  as  a  suitable  feature  for 
use  in  the  classifier  system. 

This  research  uses  a  duty  cycle  threshold  of  0.42  to  arbitrarily  decide  if  a  signal  is 
a  pulse  or  not  in  the  observed  time  interval.  The  duty  cycle  threshold  of  0.42  provided  a 
good  When  the  estimated  duty  cycle  is  less  than  or  equal  to  0.42  it  is  treated  as  a  pulse  and 
all  samples  of  \r[n] \2Smoothed  <  0.42  P0  are  set  to  zero.  This  process  allows  the  system  to 
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Figure  3.6:  Estimated  6C  over  a  range  of  SNR,/B  with  95%  confidence  intervals  calculated 
from  3,000  estimates  per  waveform  per  SNRj#  . 


ignore  samples  that  contain  only  noise  and  attempts  to  reduce  the  total  amount  of  noise  in 
the  received  observation  without  hindering  the  received  signal  pulse. 

3.3.2  Cyclic  Spectral  Correlation. 

The  SCF  for  the  received  waveforms  was  estimated  using  the  Frequency  Smoothing 
Algorithm  explained  in  Section  2.4.2.  For  brevity,  only  BPSK  and  QPSK  examples  are 
shown  to  highlight  the  cyclic  features  used  in  the  classifier.  Figure  3.7  shows  the  estimated 
SCF  for  a  simulated  BPSK  waveform  at  a  SNR  of  20  dB  and  Figure  3.8  shows  the  estimated 
SCF  for  a  simulated  QPSK  waveform  at  a  SNR  of  20  dB.  Comparing  the  figures,  one 
can  see  that  BPSK  has  a  large  cyclic  feature  for  frequency  /  =  0  and  cyclic  frequency 
a  =  2  fc  =  0.6  fs,  but  QPSK  does  not.  It  can  also  be  seen  that  the  ratio  between  the  SCF 
values  for  a  =  2  fc,  f  =  0  and  a  =  0,  /  =  fc  is  about  one  for  BPSK  and  very  low  for  QPSK. 
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(a)  Simulated  BPSK  SCF  Estimate 
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(b)  Frequency  Profiles  of  BPSK  SCF 


(c)  Cyclic  Frequency  Profiles  of  BPSK  SCF 


Figure  3.7:  Estimated  BPSK  SCF  at  SNR  =  20 dB  with  carrier  frequency,  fc  =  0.3  fS9  and  a 
bandwidth,  W  =  0.2  fs  =  2  jr— ,  using  frequency  smoothing  with  N  =  4096  and  N'  =  328. 
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(a)  Simulated  QPSK  SCF  Estimate 


(b)  Frequency  Profiles  of  QPSK  SCF 


frequency  (f) 


(c)  Cyclic  Frequency  Profiles  of  QPSK  SCF 


Figure  3.8:  Estimated  QPSK  SCF  at  SNR  =  20 dB  with  carrier  frequency,  fc  =  0.3  fs,  and  a 
bandwidth,  W  =  0.2  fs  =  2  ,  using  frequency  smoothing  with  N  =  4096  and  N'  =  328. 
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(a)  Simulated  BPSK  SCF  Estimate 


(b)  Frequency  Profiles  of  BPSK  SCF 


frequency  (f) 


(c)  Cyclic  Frequency  Profiles  of  BPSK  SCF 


(d)  Cyclic  Frequency  Profiles  of  BPSK  SCF  with  Adjusted  PSD 


Figure  3.9:  Estimated  BPSK  SCF  at  SNR  =  -5 dB  with  carrier  frequency,  fc  =  0.3  fs,  and  a 
bandwidth,  W  =  0.2  fs  =  2  ,  using  frequency  smoothing  with  N  =  4096  and  N'  =  328. 


48 


Table  3.1:  SCF  Classifier  Features 


Furthermore,  Figure  3.9  shows  the  estimated  SCF  for  a  simulated  BPSK  with  SNR  = 
-5 dB.  It  is  apparent  by  comparing  Figure  3.7  and  Figure  3.9  that  the  noise  floor  for  the 
PSD  (SCF  where  a  =  0  )  has  increased  and  the  cyclic  feature  S  (ft,  0)  has  decreased. 
To  maintain  a  ratio  of  about  1  between  the  specified  SCF  estimates,  the  noise  floor  is 
subtracted  from  the  PSD  estimate.  This  adjusted  PSD  estimate  is  denoted  as  Sx  (n,  fo)Ad] 
and  the  result  is  shown  in  Figure  3.9d.  It  is  easily  seen  that  the  ratio  between  the  S ^  (ft,  0) 
estimate  and  the  adjusted  Sx  (ft,  fc)Adj  estimate  is  about  1. 

S%  ( n,fc)Adj  =  S°Xt  ( njc )  -  Noise  Floor  (3.16) 

Four  variants  of  these  ratios  are  used  as  features  in  the  classifier  system.  One  feature, 

/  S2xfc(n,  0)  \ 

which  has  been  discussed,  is  the  ratio  between  specific  SCF  estimates  I  Sx  Inf)^-  ]•  Another 
feature  is  the  ratio  between  the  estimated  SCF  frequency  slices’  maximum  peaks  for 

(  MAXf(s2xfc(n,f ))  ) 

a  =  2 fc  and  a  =  0,  .  These  two  features  are  calculated  using  the  frequency 

smoothing  SCF  estimation  technique  from  Equation  (2.49).  Each  signal  observation  is 
N  =  4096  samples  long.  A  less  reliable,  ‘high’  resolution  SCF  estimate  where  N'  =  36 
and  a  more  reliable,  ‘low’  resolution  SCF  estimate  where  N'  =  328  are  used  to  provide 
the  four  features  used  in  the  classifier  system.  The  features  used  from  SCF  analysis  are 
explicitly  shown  in  Table  3.1.  All  the  simulated  waveforms’  values  for  the  ‘low  resolution’ 
feature  in  column  1  of  Table  3.1  over  the  SNR  range  tested  are  shown  in  Figure  3.10  with 
their  95%  confidence  intervals  for  3,000  simulated  observations  per  waveform  per  SNR. 
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Figure  3.10:  Estimated  SCF  feature  ratio  using  ‘low’  resolution  feature  from  Table  3.1 
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For  computational  efficiency,  only  the  SCF  slices  for  the  relevant  cyclic  frequencies,  a  =  0, 
and  a  =  2 fc,  are  estimated  to  calculate  the  features. 

3.3.3  Cyclic  Cumulants. 

The  received  waveform  r  [n]  cumulants  are  calculated  from  the  discrete,  time-sampled 
received  waveforms’  estimated  moments.  The  discrete  TMF  for  zero  time-lag  is  modified 
from  Equation  (2.51)  to  the  form 

Rr  [ m ,  r  =  =  E  [r  [m]n~q  (r*  [m])q]  «  r  [m]n~q  (r*  [m])q  (3.17) 

where  m  is  the  discrete  sample  index,  and  R,  [m,  r  =  ()}n  q  is  the  nth-order,  ^-conjugate 
moments  for  time  sample  m.  The  moments  were  estimated  according  to  Equation  (3.17) 
and  the  cumulant  equations  from  Table  2.5  were  used  to  calculate  the  estimated  nth-order 
^-conjugate  cumulants.  Taking  the  cumulant’s  discrete  fourier  transform  (DFT)  gives  the 
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(a)  Estimated  BPSK  2nd-Order  Cumulant  Spectrums 
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(b)  Estimated  BPSK  4th-Order  Cumulant  Spectrums 
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(c)  Estimated  BPSK  6th-Order  Cumulant  Spectrums 


Figure  3.11:  Estimated  ft-order  ^-conjugate  cyclic  cumulant  spectrums  calculated  from  a 
BPSK  waveform  with  SNR=  20 dB  using  Equation  (3.18). 


CTCF  which  is  used  to  calculated  the  cyclic  temporal  cumulants  (CTCs). 

N- 1 

Cf  [m,  0]„,9  =  £  C,  [m\  0]n  q  e~^T‘  (3.18a) 

m'= 0 

N- 1 

=  J]cx[m',  0\me~^  (3.18b) 

m'= 0 

The  estimated  CTCF  for  a  BPSK  waveform  with  SNR  =  20 dB  over  the  cyclic 
frequencies  0  <  yS  <  1  is  shown  in  Figure  3.11.  The  cyclic  cumulant’s  maximum  occurs 
at  the  cyclic  frequency  (3  =  (n-  2 q)  fc.  Also,  one  can  notice  that  secondary  peaks  occur  at 
f3  =  (n-  2 q)  fc  ±  jJ—  corresponding  to  ±  the  symbol  rate  from  the  max  peak.  Since  not  all 
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Table  3.2:  Cyclic  Cumulant  Features 


modulation  types  have  peaks  in  every  n-order  ^-conjugate  cyclic  cumulant  spectrum,  these 
cumulants  are  utilized  for  features  in  the  proposed  classifier  system.  The  peaks  of  these 
CTCFs  may  also  be  used  to  estimate  the  carrier  frequency  and  symbol  rates  of  unknown 
received  waveforms. 

For  computational  efficiency,  instead  of  calculating  the  entire  cyclic  cumulant 
spectrum,  the  classifier  system  only  estimates  the  cyclic  cumulants’  spectrums  at  the 
single  frequency  of  interest,  yS  =  in  -  2 q)  fc,  where  the  CTCF  attains  its  maximum 
value.  The  ratios  between  the  estimated  cyclic  cumulants’  magnitude  at  cyclic  frequency 
[3  =  (n  -  2 q)  fc  are  used  for  features.  The  features  calculated  from  the  estimated  CTCFs  are 
listed  in  Table  3.2. 

3.4  Classifier  Training 

After  the  received  waveform’s  features  are  estimated,  they  are  formatted  into  the 
feature  vector  if/. 

,^25]  (3.19) 

Feature  vectors  estimated  from  received  waveforms  whose  modulations  are  known  a  priori 
are  then  used  to  ‘train’  the  classifier.  This  form  of  classifier  training  is  known  as  supervised 
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Table  3.3:  Classifier  Features 


learning  and  is  shown  in  Figure  3.12.  The  classifier  uses  ‘training’  if/  vectors  for  each 
modulation  type  to  estimate  the  25-dimensional  PDF  for  that  modulation.  All  features  used 
by  the  classifier  system  are  listed  in  Table  3.3. 

The  classification  model  used  in  this  research  assumes  that  the  <A  feature  values  for 
each  class  have  a  multivariate-Gaussian  distribution  represented  by 

p  (W)  = - - TexP  (*A  -  Pk)  £ k 1  Uf  ~  fik)  (3.20) 

V-  ;  (27r)T  det(^)1  2  V“  V“ 
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Figure  3.12:  Classifier  Training 

where  P  (ij/\kj  is  the  PDF  of  a  feature  vector,  ip,  for  a  given  signal  class  k,  E*  is  the  estimated 
feature’s  sample  covariance  matrix  for  signal  k,  juk  is  a  vector  of  estimated  feature  sample 
means  of  signal  k,  (-)T  denotes  the  matrix  transpose  operation,  and  E"1  is  the  inverse  matrix 
ofE*. 

The  sample  mean  and  unbiased  sample  covariance  for  signal  k  are  estimated  during 
supervised  training  from  M  independent  estimations  of  signal  /c’s  features 

(3.21a) 

m=  1 

1  M  T 

%  =  M  _  {  At)  (*h  ~  fa)  (3.21b) 

m- 1 

where  ipk  are  calculated  features  from  observations  known  to  be  from  signal  class  k. 

The  classifier  was  trained  for  each  signal  class  at  each  SNR  using  <A  vectors  extracted 
from  1,000  independently  simulated  waveforms  of  each  signal  class  at  each  simulated 
SNR.  All  twelve  signal  class  waveforms  were  simulated  for  26  levels  of  SNR  from 
20 dB  to  -5dB,  totaling  12  x  1,000  =  12,000  tfj  vectors  used  for  training  each  SNR  and 
12,000  x  26  =  312,000  \J/  vectors  used  for  training  overall.  The  MATLAB®  command 
‘ClassificationDiscriminant.fit’  with  the  ‘diagQuadratic’  option  was  used  to  estimate  the 
signal  class  feature  PDFs,  P  (ip\k}  from  the  training  observations. 
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Figure  3.13:  Test  the  Classifier 


3.5  Performance  Criteria 


Classifier  predictions  from  independent  test  observations  are  used  to  approximate  the 
classifier  performance  for  each  simulated  SNR.  Using  the  estimated  P  (ifr\ kj  for  each  signal 
class  developed  through  the  training  process  shown  in  Figure  3.12,  the  classifier  model 
calculates  the  estimated  posterior  probabilities 


P(k\l) 


P(<iAk)p(k)  p(ijAk)P(k) 


(3.22) 


2  PU\k')P(k’) 

k'= 1  V—  ' 

where  P  (k\if/j  is  the  probability  that  an  estimated  if/  vector  is  from  signal  k,  the  denominator 
term  is  a  normalization  constant  to  assure  that  2  P(k'\<p)  =  1,  and  P  (k)  is  the  prior 
probability  of  if/  belonging  to  signal  class  k.  For  this  classifier,  all  signal  classes  are  assumed 
to  be  equally  likely  so  P(k)  =  j  for  all  signal  classes.  This  is  a  classifier  predicts  the  signal 
class  that  if/  belongs  to  based  on  the  maximum  a  posteriori  (MAP)  likelihood 


k  =  arg^max  p(k\if/j  (3.23) 

The  classifier’s  performance  is  assessed  through  the  process  shown  in  Figure  3.13  with 
3,000  independently  simulated  observations  of  each  signal  class  at  each  known  SNR.  The 
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Figure  3.14:  Confusion  Matrix 


25  features  are  estimated  from  each  observation  and  formatted  into  the  feature  vector  if/, 
if/  is  then  the  input  to  the  trained  classifier  which  calculates  posterior  probabilities  that 
the  estimated  if/  belongs  to  the  signal  class  k.  The  posterior  probability  for  the  signal 
class  k  with  the  maximum  likelihood  is  then  selected  as  the  predicted  signal  class  k.  The 
MATLAB®  ‘predict’  command  is  used  to  calculate  and  compare  the  posterior  probabilities 
(P  to  make  the  MAP  likelihood  classification  decision  in  Equation  (3.23). 


The  ratio  of  correct  class  predictions  over  total  predictions  is  the  probability  of  correct 
classification 


%C  = 


Total  #  of  trials  where  kT  =  k 


(3.24) 


Total  #  of  trials 

where  kT  represents  the  true  signal  class  and  k  is  the  predicted  signal  class.  This  ratio  will 
represent  the  classifier’s  overall  performance.  Confusion  matrices  and  receiver  operating 
characteristic  (ROC)  curves  will  also  be  used  to  illustrate  the  classifier’s  performance  for 
each  signal  class  at  SNRs  of  interest.  A  confusion  matrix  provides  the  true  class  and 
predicted  class  information  for  all  trials  as  shown  in  Figure  3.14.  A  ROC  curve  illustrates 
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Ideal 


the  probability  of  detection  ( PD )  versus  the  probability  of  false  alarm  (P  fa)  for  a  varying 
threshold  A  at  a  specific  SNR.  For  this  research’s  application,  PD  is  the  probability  that 
the  true  signal  class  will  be  correctly  detected  with  a  threshold  A  and  PFA  is  the  probability 
that  an  incorrect  signal  class  will  be  falsely  detected  as  the  true  signal  class  for  the  same 
A.  These  probabilities  are  calculated  from  the  posterior  probability  mass  functions  (PMFs) 
generated  by  classifier  test  results  using 

Total  #  of  trials  where  [p  (k  =  kT\ip^  >  d] 


Pd,  a  - 


Total  #  of  trials  where  k  =  kT 


FA,  A 


Total  #  of  trials  where  [  p  (k  ±  kT\i//j  >  A 
Total  #  of  trials  k±kT 


(3.25a) 


(3.25b) 


Examples  of  generic  ROC  curves  are  shown  in  Figure  3.15. 
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IV.  Results  and  Analysis 


This  chapter  will  describe  the  simulations,  present  results,  and  analyze  the  modulation 
classification  system’s  performance.  It  is  done  by  generating  probability  of  correct 
classification  versus  SNR^g  graphs,  ROC  curves,  and  confusion  matrices.  The  feature- 
based  modulation  classification  system’s  ability  to  recognize  various  signal  class  groupings 
will  also  be  tested. 

In  Section  4.1  the  test  simulations  are  explained,  including  all  pertinent  signal 
parameters  used  for  both  training  and  testing.  Section  4.2  presents  the  test  simulation 
results  and  provides  an  analysis  of  them. 

4.1  Simulation  Setup 

The  designed  modulation  classification  system  was  tested  for  the  signal  modulations: 
BPSK,  QPSK,  16-QAM,  64-QAM,  8-PSK,  16-PSK,  Bi-phase  Barker5,  Bi-phase  Barkern, 
Bi-phase  Barkers  n,  Frank49,  PX49,  and  LFM.  First,  each  signal  modulation  was 
independently  simulated  4,000  times  for  each  of  26  SNRrfg  levels,  a  pulse  shaping  filter  with 
50%  excess  bandwidth  was  used  to  band-limit  the  signals,  and  independent  realizations  of 
AWGN  were  added  to  each  simulated  signal  modulation  to  achieve  the  desired  simulated 
SNRjg.  Each  simulated  waveform  is  212  =  4,096  samples  long  with  a  normalized 
intermediate  frequency  (IF)  bandwidth  of  W  =  2B  =  0.2 fs  and  was  modulated  to  a  carrier 
frequency  fc  =  0.3  fs.  It  is  assumed  that  the  carrier  frequency  and  the  SNR  have  been 
accurately  determined  and  are  available  to  the  modulation  classification  system  for  feature 
extraction  and  classification. 

The  first  test  simulation  provides  a  baseline  for  the  modulation  classification  system’s 
performance.  25  signal  features  were  estimated  from  1,000  observations  of  each  signal 
modulation  at  all  SNRrfs  levels  and  were  used  for  training  the  classifier.  This  process 
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was  described  in  Section  3.4  and  shown  in  Figure  3.12.  Then,  the  25  features  were 
estimated  from  3,000  new  observations  of  each  signal  modulation  and  were  used  to  test 
the  classifier’s  performance  at  each  simulated  SNR(/B.  The  testing  process  was  described  in 
Section  3.5  and  illustrated  in  Figure  3.13.  The  classifier’s  performance  for ‘ideal’  conditions 
is  evaluated  by  using  training  and  testing  observations  with  the  same  critical  parameters 
such  as  bandwidth  and  carrier  frequency. 

This  ideal  case  is  assessed  for  three  classification  scenarios.  First,  the  scenario  where 
the  system  must  classify  all  distinct  signal  modulations  that  were  simulated.  Second, 
the  signal  modulations  are  grouped  by  their  general  modulation  family  type  to  assess  the 
system’s  ability  to  recognize  and  classify  a  signal’s  modulation  family.  BPSK,  QPSK,  and 
LFM  remain  their  own  distinct  classes,  but  16-QAM  and  64-QAM  are  labeled  collectively 
as  M-QAM;  8-PSK  and  16-PSK  are  labeled  collectively  as  M-PSK;  Barker5,  Barkern, 
and  Barker5  i  |  are  labeled  collectively  as  Bi-Phase  Barker;  and  last,  Frank49  and  Px49  are 
labeled  collectively  as  Poly  Phase.  Third,  the  signal  modulations  are  grouped  into  the  broad 
categories  of  communication  and  pulsed  radar  waveforms.  The  communication  category 
includes  BPSK,  QPSK,  16-QAM,  64-QAM,  8-PSK,  and  16-PSK  while  the  pulsed  radar 
category  includes  Bi-phase  Barker5,  Bi-phase  Barkern,  Bi-phase  Barker5>n,  Frank49,  Px49, 
and  LFM.  The  classifier  is  re-trained  and  re-tested  for  these  relaxed  signal  groupings 
to  assess  its  performance  when  only  general  information  about  an  unknown  waveform’s 
modulation  is  required. 

Next,  the  modulation  classification  system’s  robustness  to  carrier  frequency  and 
bandwidth  deviation  between  the  observations  used  for  training  and  the  observations  used 
for  testing  were  assessed.  The  classification  system  is  still  trained  with  1,000  observations 
and  tested  with  3,000  observations  of  each  signal  modulation  at  all  SNR^g  levels.  The 
training  observations  were  simulated  with  W  =  2B  =  0.2 fs  and  fc  =  0.3  fs  for  both  tests, 
but  the  test  observations  had  varied  simulated  signal  parameters;  first  with  fc  =  0.3  fs  and 
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an  increased  IF  bandwidth  W  =  2B  =  OAfs,  and  then  with  IF  bandwidth  W  =  2B  =  0.2 fs 
and  a  decreased  carrier  frequency  fc  =  0.2 fs . 

Probability  of  correct  classification  (%C)  versus  SNR,/B  graphs,  ROC  curves,  and 
confusion  matrices  are  generated  for  each  simulation  test.  These  measures  of  performance 
were  explained  in  detail  in  Section  3.5  and  are  used  to  attain  a  comprehensive  assessment 
for  the  modulation  classification  system  developed  in  this  research  work. 

4.2  Classifier  Performance  with  Ideal  Training  Data 

4.2.1  Signal  Modulation  Type  Classification. 

The  designed  modulation  classification  system  was  first  simulated  for  a  classification 
scenario  with  ideal  conditions.  Figure  4.1  shows  the  performance  versus  simulated  SNR^g 
for  the  %C  modulation  classification  for  all  twelve  signal  modulations  considered:  BPSK, 
QPSK,  16-QAM,  64-QAM,  8-PSK,  16-PSK,  Bi-phase  Barker5,  Bi-phase  Barkern,  Bi¬ 
phase  Barker5  i  i ,  Franks,  Px49,  and  LFM.  This  is  broken  down  into  the  classifier’s 
classification  performance  for  each  simulated  signal  modulation  in  Figure  4.2.  Per  the 
simulation  results  in  Figure  4.1  and  Figure  4.2,  the  modulation  classification  system  has  an 
overall  average  %C  >  90%  for  the  simulated  signal  modulations  tested  at  SNR  =  9dB. 
Table  4.1  provides  the  SNR  =  9dB  confusion  matrix  for  the  system’s  performance  in 
Figure  4.2. 

Figure  4.2a  and  Figure  4.2c  shows  that  classifier  performances  for  BPSK,  QPSK, 
16-QAM,  64-QAM,  LFM,  Barker5,  Barkern,  and  Barker5>11  are  generally  increasing  over 
SNR.  However,  there  is  some  confusion  in  the  simulated  classifier  system  for  the  other 
signal  modulation  types  at  low  SNR  levels.  The  performance  for  8-PSK  and  16-PSK  seem 
to  be  negatively  correlated  in  Figure  4.2b  up  to  SNR  =  -4dB  and  the  performance  for 
Frank49  and  Px49  seem  to  be  negatively  correlated  in  Figure  4.2d  up  to  SNR  =  8dB. 
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Table  4. 1 :  Signal  Modulation  Type  Confusion  Matrix  for  SNR  =  9dB  in  Figure  4. 1 


Communication 

Pulse  Compression  Radar 

BPSK 

QPSK 

16QAM 

64QAM 

8PSK 

16PSK 

Barkers 

Barkern 

Barker5  n 

Frank49 

Px49 

LFM 

Communication 

BPSK 

3000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

QPSK 

0 

3000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

16QAM 

0 

0 

2777 

223 

0 

0 

0 

0 

0 

0 

0 

0 

64QAM 

1 

0 

143 

2855 

0 

0 

0 

1 

0 

0 

0 

0 

8PSK 

0 

0 

0 

0 

2999 

1 

0 

0 

0 

0 

0 

0 

16PSK 

0 

0 

0 

0 

3 

2997 

0 

0 

0 

0 

0 

0 

Pulse  Compression 
Radar 

Barker5 

0 

0 

0 

0 

0 

0 

2827 

173 

0 

0 

0 

0 

Barkern 

0 

0 

0 

0 

0 

0 

391 

2607 

2 

0 

0 

0 

Barker511 

0 

0 

0 

0 

0 

0 

2 

5 

2993 

0 

0 

0 

Frank49 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2496 

499 

5 

Px49 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1771 

1214 

15 

LFM 

0 

0 

0 

0 

1 

0 

0 

0 

0 

2 

5 

2992 

SNRdB 


Figure  4.1:  Classifier  average  performance,  for  all  12  signal  modulation  types  considered, 
developed  from  48,000  simulated  observations  with  W  =  0.2 fs  and  fc  =  0.3 fs.  12,000 
observations  were  used  for  training,  1,000  from  each  signal  modulation  type,  and  36,000 
observations  were  used  for  testing,  3,000  from  each  signal  modulation  type. 


61 


(a)  BPSK,  QPSK,  16-QAM,  64-QAM,  and  LFM  (b)  8-PSK  and  16-PSK 


Figure  4.2:  Specific  modulation  type  classification  performance  of  system  in  Figure  4.1 

At  SNR  =  9dB,  the  simulated  system’s  overall  average  classification  performance 
%C  is  greater  than  90%,  but  Table  4.1  shows  that  more  than  half  of  the  PX49  modulated 
pulses  were  misclassified  as  Frank49  modulated  pulses  and  about  16%  of  the  Franl<49 
modulated  pulses  were  misclassified  as  Px49  modulated  pulses.  This  should  not  be 
too  surprising  since  Franl^  and  PX49  coded  modulated  pulses  have  very  similar  phase 
sequences.  Additionally,  all  other  simulated  signal  modulation  types,  besides  Barker!  1, 
boast  less  than  10%  misclassification  and  BPSK,  QPSK,  8-PSK,  16-PSK,  Barker5>n,  and 
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Figure  4.3:  Classifier  System  ROCs  for  the  12  Modulation  Types  at  SNR  =  9dB 

LFM  were  classified  with  greater  than  99%  accuracy.  It  is  also  interesting  to  note  that  only 
one  simulated  observation  of  pulse  compression  radar  modulation  was  misclassified  as  a 
communication  modulation  and  vice  versa.  Also,  the  misclassifications  are  mainly  between 
different  orders  of  similar  modulation  types  such  as  16-QAM  and  64-QAM;  8-PSK  and  16- 
PSK;  Barkers,  Barkern,  and  Barker5;n;  and  Franl^  and  Px49. 

Figure  4.3  shows  the  classification  system’s  ROC  curves  for  all  signal  modulations 
at  SNR  =  9dB.  The  figure  region  for  PFA  from  0  to  0.15  and  PD  from  0.85  to  1 
has  been  expanded  to  ease  figure  readability.  Remember,  that  the  ROC  curve  for  any 
signal  modulation  ‘X’,  represents  the  probability  of  signal  modulation  ‘X’  being  correctly 
detected  ( PD )  and  the  corresponding  probability  any  other  modulation  is  incorrectly 
detected  as  modulation  ‘X’  ( PFA )  for  a  sliding  threshold.  Therefore,  it  can  be  seen  in 
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Table  4.2:  Signal  Modulation  Type  Confusion  Matrix  for  SNR  =  OdB  in  Figure  4.1 


Communication 

Pulse  Compression  Radar 

BPSK 

QPSK 

16QAM 

64QAM 

8PSK 

16PSK 

Barker5 

Barkern 

Barker511 

Frank49 

Px49 

LFM 

Communication 

BPSK 

2991 

0 

0 

0 

0 

0 

0 

9 

0 

0 

0 

0 

QPSK 

0 

2914 

79 

1 

6 

0 

0 

0 

0 

0 

0 

0 

16QAM 

0 

71 

2147 

731 

28 

14 

0 

5 

0 

3 

1 

0 

64QAM 

0 

14 

1509 

1389 

24 

40 

0 

15 

0 

5 

4 

0 

8PSK 

0 

0 

48 

33 

2109 

802 

7 

1 

0 

0 

0 

0 

16PSK 

0 

1 

53 

33 

2089 

816 

7 

1 

0 

0 

0 

0 

Pulse  Compression 
Radar 

Barker5 

0 

0 

0 

0 

0 

0 

808 

1696 

182 

54 

228 

32 

Barkern 

0 

0 

0 

0 

0 

0 

187 

1960 

526 

78 

193 

56 

Barker511 

0 

0 

0 

0 

0 

0 

15 

219 

2762 

0 

4 

0 

Frank49 

0 

0 

0 

0 

0 

0 

8 

25 

2 

712 

2152 

101 

x" 

CL 

0 

0 

0 

0 

0 

0 

11 

26 

3 

643 

2207 

110 

LFM 

0 

0 

0 

0 

0 

0 

24 

60 
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1 
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Figure  4.4:  Classifier  System  ROCs  for  the  12  Modulation  Types  at  SNR  =  OdB 
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Figure  4.3  that  specific  thresholds  can  be  set  to  achieve  PD  =  90%  with  PFA  <  2% 
for  every  modulation  type  but  PX49  and  Frank49  in  a  9dB  SNR  environment  with  this 
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simulated  modulation  classification  system.  The  simulated  PX49  and  Frank49  modulations 
must  sacrifice  PFA  ~  7%  to  achieve  PD  =  90%. 

Even  though  the  simulated  classification  system  achieved  %C  <  60%  in  a  OdB 
SNR  environment,  Table  4.2  and  Figure  4.2a  shows  that  BPSK  and  QPSK  can  still  be 
reliably  recognized  by  the  classification  system  with  %C  >  95%  accuracy  in  such  noisy 
environments.  Also,  per  Figure  4.4,  the  classification  system  can  achieve  PD  =  90%  with 
PFA  <  10%  for  the  simulated  communication  modulations  and  PD  =  90%  with  PFA  <  20% 
for  the  simulated  pulse  compression  radar  modulations  by  using  specific  thresholds. 

4.2.2  Signal  Modulation  Family  Classification. 

For  this  simulation  test,  the  system’s  performance  recognizing  and  classifying  general 
modulation  families  is  analyzed.  Using  the  same  simulated  48,000  observations  at  every 
SNR  as  before,  the  classifier  was  trained  and  tested  to  recognize  the  modulation  families 
given  in  Table  4.3.  Figure  4.5  shows  this  simulated  classifier’s  overall  average  probability 
of  correct  classification  (%C)  versus  SNR^g  performance  for  the  seven  modulation  families 
considered  compared  to  the  previous  classifier  for  modulation  type  classification.  The 


Table  4.3:  Modulation  Families 


Modulation  Family 

Modulation  Types 

#  of  Training/Test 
Observations  per  SNR 

BPSK 

BPSK 

1,000/3,000 

QPSK 

QPSK 

1,000/3,000 

M-QAM 

16-QAM,  64-QAM 

2,000/6,000 

M-PSK 

8-PSK,  16-PSK 

2,000/6,000 

Bi-Phase  Barker 

Barker5,  Barkern, 
Barker511 

3,000/9,000 

Poly  Phase 

Frank49,  Px49 

2,000/6,000 

LFM 

LFM 

1,000/3,000 
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Figure  4.5:  Classifier  average  performance,  for  the  7  modulation  families  given  in 
Table  4.3,  developed  from  48,000  simulated  observations  with  W  =  0.2 fs  and  fc  =  0.3  fs. 
12,000  simulated  observations  were  used  for  training  and  36,000  simulated  observations 
were  used  for  testing. 


classifier’s  classification  performance  for  each  simulated  modulation  family  is  shown  in 
Figure  4.6. 

Per  the  simulation  results  in  Figure  4.5  and  Figure  4.6,  this  modulation  classification 
system  achieves  overall  %C  «  88%  at  SNR  =  OdB  and  the  simulated  modulation  families 
can  be  correctly  identified  with  overall  %C  =  90%  at  an  SNR  about  8dB  less  than  that 
required  for  simulated  modulation  type  classification.  Also,  every  modulation  family  has  a 
%C  >  90%  for  SNR  >  6dB. 

Table  4.4  provides  the  SNR  =  OdB  confusion  matrix  for  the  system  in  Figure  4.5. 
Since  some  general  modulation  families  encompass  multiple  specific  modulation  types, 
the  total  number  of  observations  used  for  training  and  testing  each  modulation  family  for 
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Figure  4.6:  Modulation  family  classification  performance  of  system  in  Figure  4.5  with 
overall  average  system  classification  performance 

this  simulation  is  also  given  in  Table  4.3.  Nearly  85%  of  LFM  modulated  pulses  were 
misclassified  as  Poly  Phase  modulated  pulses,  but  the  signal  modulations  in  the  other 
modulation  families  are  recognized  with  %C  >  89%.  Also,  there  was  only  one  simulated 
pulse  compression  radar  modulation  family  observation  misclassified  as  a  communication 
modulation  family  and  only  51  of  18,000  simulated  communication  modulation  family 
observations  were  misclassified  as  a  pulse  compression  radar  modulation  family. 

The  ROC  curves  for  the  modulation  family  classification  system  at  SNR  =  OdB  in 
Figure  4.5  are  given  in  Figure  4.7.  The  modulation  family  classification  system  can  set 
thresholds  to  achieve  PD  =  90%  and  P fa  <  1%  for  all  but  the  simulated  Poly  Phase  and 
LFM  signals  at  this  SNR.  The  simulated  Poly  Phase  pulse  compression  modulations  can 
achieve  PD  =  90%  and  PFA  ~  10%  and  the  simulated  LFM  can  achieve  PD  =  90%  and 
P fa  *  17%. 
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Figure  4.7:  Classifier  System  ROCs  for  7  Modulation  Families  at  SNR  =  OdB 


Table  4.4:  Modulation  Family  Confusion  Matrix  for  SNR  =  OdB  in  Figure  4.5 


Communication 

Pulse  Compression  Radar 

#Test 

Observations 

BPSK 

QPSK 

MQAM 

MPSK 

Bi-Phase 

(Barker) 

Poly  Phase 
(Frank,  Px) 

LFM 

Communication 

BPSK 

2991 

0 

0 

0 

9 

0 

0 

3,000 

QPSK 

0 

2920 

75 

5 

0 

0 

0 

3,000 

MQAM 

0 

100 

5738 

124 

12 

26 

0 

6,000 

MPSK 

0 

1 

167 

5828 

4 

0 

0 

6,000 

Pulse  Compression 
Radar 

Bi-Phase 

(Barker) 

1 

0 

0 

0 

8075 

741 

183 

9,000 

Poly  Phase 
(Frank,  Px) 

0 

0 

0 

0 

21 

5675 

304 

6,000 

LFM 

0 

0 

0 

0 

46 

2500 

454 

3,000 

68 


o 


c 

o 

c5 

g 

C/D 

c/) 

_ro 

O 

o 

0 

o 

O 

> 


■Q 

03 

X) 

O 


1 

0.98 
0.96 
0.94 
0.92 
0.9 

0.88^ 
0.86 
0.84 
0.82 


0.8 


© 


-5 


0 


0' 


© 


6 

i 


.o' 


© 


p0 


0 


■0 


-©- 


Radar  vs  Comm 
Classification 


q  _  Modulation  Family 
Classification 

q  _  Modulation  Type 
Classification 


10 


15 


SNR 


dB 


Figure  4.8:  Classifier  average  performance,  for  distinguishing  between  communication 
and  pulsed  radar  modulations  according  to  Table  4.5.  System  is  developed  from  48,000 
simulated  observations  with  W  =  0.2 fs  and  fc  =  0.3 fs.  A  total  of  12,000  simulated 
observations  were  used  for  training  and  36,000  simulated  observations  were  used  for 
testing. 


4.2.3  Communication  vs.  Pulse  Compression  Radar  Modulation  Classification. 

Last,  the  system’s  ability  to  distinguish  the  simulated  communication  modulations 
from  the  simulated  radar  pulse  compression  modulations  is  shown.  For  this  simulation, 
the  classification  system  is  trained  and  tested  for  classifying  the  12  signal  modulations  as 
either  a  communication  modulation  or  pulse  compression  radar  modulation  as  illustrated 
in  Table  4.5.  Again,  the  same  48,000  observations  at  every  SNR^s  (4,000  for  each  of 
12  signal  modulation  types)  were  used  to  train  and  test  the  classification  system  for 
this  binary  classification  decision.  Figure  4.8  shows  the  simulated  pulse  radar  versus 
communication  classifier’s  overall  probability  of  correct  classification  (%C)  versus  SNR,/# 
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Table  4.5:  Radar  and  Communication  Waveforms 


Signal  Type 

Modulation  Types 

#  of  Training/Test 
Observations  per  SNR 

Communication 

BPSK,  QPSK,  16-QAM, 
64- QAM,  8-PSK,  16-PSK 

6,000/18,000 

Barker5,  Barkern/ 

Radar 

Barker5  n,  Frank49/ 

6,000/18,000 

Px49/  LFM 

Figure  4.9:  Pulsed  radar  and  communication  classification  performance  of  system  in 
Figure  4.8  with  overall  average  system  classification  performance 


performance  compared  to  the  previous  classifiers  simulated,  and  Figure  4.9  shows  the 
system’s  individual  %C  for  pulse  compression  radar  and  communication  modulations. 

Per  the  simulation  results  in  Figure  4.8,  this  modulation  classification  system  achieves 
%C  >  98%  at  SNR  =  -ldB  and  even  maintains  %C  >  90%  at  SNR  =  -4dB.  It 
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Figure  4.10:  Classifier  system  ROCs  for  communication  detection  at  SNR  =  -5 dB.  The 
ROC  for  pulsed  radar  detection  is  simply  the  reflection  of  the  communication  detection 
ROC  about  Pd  -  1  -  Pfa- 

can  achieve  overall  %C  =  90%  between  pulse  compression  radar  and  communication 

modulations  at  an  SNR  5dB  lower  than  it  can  for  modulation  family  classification  and 

about  12dB  lower  than  modulation  type  classification.  The  simulations  show  in  Figure  4.8 

that  the  proposed  system  can  provide  reliable  %C  >  90%  for  modulation  types  down  to 

SNR  «  8dB,  modulation  families  down  to  SNR  «  0.5dB,  and  pulse  compression  radar 

versus  communication  modulation  down  to  SNR  «  -4dB. 

In  Figure  4.9,  it  is  seen  that  pulse  compression  radar  modulations  are  classified  with 

%C  less  than  the  communication  modulation  waveforms  and  never  achieve  %C  >  99%. 

Therefore,  the  overall  %C  for  the  system  never  achieves  %C  =  100%.  The  communication 

modulations  are  classified  with  very  good  reliability  maintaining  %C  >  96%  for  SNR  > 

-5dB,  but  the  pulse  compression  radar  modulation  %C  degrades  to  %C  ~  78%  by 

SNR  =  -5dB. 
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Table  4.6:  Communication  vs  Pulsed  Radar  Modulation  Confusion  Matrix  for  SNR  = 
-ldB  in  Figure  4.8 


Communication 

Pulse  Compression 
Radar 

Communication 

17914 

86 

Pulse  Compression 
Radar 

226 

17774 

Table  4.7:  Communication  vs  Pulsed  Radar  Modulation  Confusion  Matrix  for  SNR  = 
-5dB  in  Figure  4.8 


Communication 

Pulse  Compression 
Radar 

Communication 

17503 

497 

Pulse  Compression 
Radar 

3865 

14135 

Since  this  test  scenario  is  a  binary  classification  decision,  only  the  communication 
modulation  ROC  curve  is  shown  in  Figure  4.10.  The  pulsed  radar  modulation  ROC  curve 
can  be  inferred  from  Figure  4.10  as  the  reflection  about  PD  =  1  -  PFA.  Thresholds  in 
this  simulated  classifier  system  can  be  set  to  achieve  PD  >  90%  for  both  communication 
and  pulsed  radar  modulations  with  PFA  <  10%  at  SNR  =  -5dB.  The  confusion  matrix 
results  for  the  system  in  Figure  4.8  at  SNR  =  -ldB  and  SNR  =  -5dB  are  given  in 
Table  4.6  and  Table  4.7  respectively.  From  both  Table  4.6  and  Table  4.7,  one  can  see  that 
the  system  tends  to  incorrectly  classify  pulse  compression  radar  modulations  more  than 
communication  modulations. 
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Figure  4.11:  Classifier  performance  with  mismatched  bandwidth  between  training  and  test 
observations  for  all  12  signal  modulations  considered.  Trained  from  12,000  simulated 
observations,  1,000  for  each  signal  modulation  type,  with  W  =  0.2/*  and  fc  =  0.3  /  and 
tested  with  36,000  simulated  observations,  3,000  from  each  signal  modulation  type,  with 
fc  =  0.3/*  and  an  increased  W  =  0.4/*. 


4.3  Classifier  Bandwidth  Sensitivity 

The  modulation  classification  system’s  sensitivity  to  bandwidth  deviation  between 
the  observations  used  for  training  and  the  observations  used  for  testing  is  assessed  in  this 
section.  1 ,000  observations  of  each  signal  modulation  at  all  SNR  levels  with  W  =  2B  = 
0.2/*  and  fc  =  0.3/*  are  used  for  training,  and  the  classifier  is  tested  with  3,000  signal 
modulation  observations  at  each  SNR  with  the  same  fc  =  0.3/* ,  but  an  increased  bandwidth 
W  =  2B  =  0.4/*.  The  modulation  classification  system’s  sensitivity  to  mismatched 
bandwidth  between  training  and  testing  signal  observations  is  shown  by  Figure  4.11, 
Figure  4.12,  Table  4.8,  and  Figure  4.13. 
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(a)  BPSK,  QPSK,  16-QAM,  64-QAM,  and  LFM 


(c)  Barkers,  Barkern,  and  Barkers^ 


Figure  4.12:  Specific  signal  modulation  classification  performances  compared  between  the 
systems  with  overall  average  classification  performance  shown  in  Figure  4.1 1.  Subscripted 
‘BW’  indicates  mismatched  bandwidths  between  training  data  and  testing  observations 
and  the  square  markers  with  dashed  lines  are  the  corresponding  %C  performance  for 
mismatched  bandwidth. 


From  Figure  4.11  it  is  seen  that  bandwidth  has  some  effect  as  overall  classification 
performance  was  degraded  when  the  test  observations  had  different  information  bandwidths 
than  the  training  observations.  This  mismatch  in  bandwidth  had  the  largest  effect  on  the 
Frank49,  Px49,  and  Barkers  n  classification  performance.  The  confusion  matrix  for  this 
simulation  at  SNR  =  8dB  is  shown  in  Table  4.8  and  the  ROC  is  shown  in  Figure  4.13. 
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Table  4.8:  Confusion  matrix  for  the  12  modulation  types  using  mismatched  training  and 
testing  bandwidths  at  SNR  =  8dB  in  Figure  4.11 
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Figure  4.13:  Classifier  System  ROCs  for  the  12  Modulation  Types  at  SNR  =  8dB  with 
Mismatched  Bandwidth  of  W  =  OAfs 


Even  though  the  %C  using  MAP  likelihood  is  very  low  for  some  modulation  types,  the 
ROCs  in  Figure  4.13  suggests  there  is  a  threshold  for  each  modulation  type  that  will  allow 
PD  >  90%  with  PFA  <  10%  at  SNR  =  8dB  when  all  modulations  are  equally  likely. 
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Figure  4.14:  Classifier  performance  with  mismatched  carrier  frequency  between  training 
and  test  observations  for  all  12  signal  modulations  considered.  Trained  from  12,000 
simulated  observations,  1,000  for  each  signal  modulation  type,  with  W  =  0.2  fs  and  fc  = 
0.3 fs  and  tested  with  36,000  simulated  observations,  3,000  from  each  signal  modulation 
type,  with  matched  W  =  0.2 fs  and  a  decreased  fc  =  0.2 fs . 


4.4  Classifier  Carrier  Frequency  Sensitivity 

The  modulation  classification  system’s  sensitivity  to  carrier  frequency  deviation 
between  the  training  and  testing  observations  is  also  assessed.  1,000  observations  of  each 
signal  modulation  at  all  SNR  levels  with  W  =  2B  =  0.2 fs  and  fc  =  0.3  fs  were  used  for 
classifier  training  and  3,000  signal  modulation  observations  at  each  SNR  with  the  same 
bandwidth,  W  =  2B  =  0.2 fs,  but  a  decreased  carrier  frequency  fc  =  0.2 fs  were  used  for 
testing.  The  performance  sensitivity  to  carrier  frequency  is  summarized  by  Figure  4.14. 

From  Figure  4.14  it  is  seen  that  carrier  frequency  has  little  effect  on  the  modulation 
feature  PDFs  because  the  overall  %C  is  the  same  when  test  observations  had  different 
carrier  frequencies  than  the  training  observations.  There  was  also  only  little  effect  on  the 
ROCs  and  confusion  matrices  so  they  are  not  presented. 
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V.  Conclusions 


The  goal  of  this  research  was  to  advance  the  application  of  modulation  classification 
by  developing  and  simulating  a  reliable  automatic  modulation  recognition  system 
capable  of  discerning  between  a  wide  range  of  both  common  pulse  compression  radar  and 
communication  modulations.  The  performance  was  simulated  for  varying  classification 
groupings,  signal  bandwidths,  and  signal  carrier  frequencies  for  a  wide  range  of  received 
SNR. 


5.1  Summary 

A  modulation  classification  system  based  on  estimated  duty  cycle  and  cyclic  features 
was  developed  and  its  performance  characterized  for  three  distinct  class  groupings  in 
with  ‘ideal’  training  data,  where  received  testing  signal  parameters  equal  training  signal 
parameters.  Also,  performance  when  the  received  signal’s  bandwidth  or  carrier  frequency 
was  different  from  the  training  observations’  was  assessed  to  determine  system  robustness 
when  suboptimal  training  data  is  available. 

Using  the  developed,  feature-based  maximum  a  posteriori  (MAP)  classification 
system  with  ideal  training  observations,  an  average  correct  modulation  type  classification 
rate  of  %C  >  90%  was  achieved  at  SNR  =  8dB.  Furthermore,  the  BPSK  and  QPSK 
signals  were  classified  at  %C  >  95%  accuracy  at  SNR  =  0  dB;  the  classification 
system  had  the  most  trouble  distinguishing  between  the  Frank49  and  Px49  radar  pulse 
compression  modulations  given  they  have  nearly  identical  phase  sequences.  The  system’s 
performance  was  %C  ~  90%  at  SNR  =  OdB  when  it  was  re-trained  to  recognize  modulation 
families;  it  had  the  most  trouble  distinguishing  between  the  LFM  and  Poly-Phase  radar 
modulation  families.  Last,  the  system’s  performance  was  %C  ~  90%  at  SNR  =  -4dB 
when  it  was  re-trained  to  distinguish  only  between  the  pulsed  radar  and  communication 
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modulations.  As  the  noise  corrupting  the  received  signal  becomes  dominant  (lower 
SNRs),  the  cyclic  features  used  for  classification  become  less  reliable  and  classification 
performance  degrades. 

The  classification  system’s  performance  suffers  when  the  received  signal  bandwidth 
differs  from  the  training  signal’s  bandwidth.  The  system  is  still  able  to  classify  the 
communication  modulations  with  comparable  reliability,  but  many  pulse  compression  radar 
modulation’s  classification  reliability  is  greatly  affected.  When  training  signal’s  bandwidth 
W  =  0.2  fs  and  testing  signal’s  bandwidth  W  =  0.4  fs,  the  overall  %C  was  degraded  by  about 
10%-14%  over  the  simulated  SNR  range.  Alternatively,  the  mismatched  carrier  frequency 
simulation  showed  little  effect  on  the  system’s  performance.  Provided,  accurate  carrier 
frequency  information  is  available  to  the  classification  system,  the  mismatched  carrier 
frequency  did  not  effect  the  classification  feature  PDFs.  These  simulations  indicated  that 
the  developed  classification  system  may  be  robust  to  variances  in  carrier  frequency,  but 
requires  matched  bandwidth  between  training  signal  data  and  received  signal’s  to  provide 
reliable  classification  performance. 

5.2  Impact 

This  research  has  taken  current  communication  modulation  classification  methods 
and  applied  them  to  a  larger  signal  set  that  includes  both  communication  and  pulse 
compression  radar  modulations.  It  has  also  incorporated  many  features  shown  to  be 
successful  in  current  literature,  such  as  duty  cycle,  spectral  correlation,  and  cyclic 
cumulants,  into  one  classification  system.  The  results  provide  evidence  that  an  integrated 
modulation  classification  system  with  application  for  both  communication  and  pulsed  radar 
modulations  is  feasible  for  applications  in  both  Cognitive  SDR,  Cognitive  Radar,  and  EW 
operations. 

A  very  broad  range  of  applications  have  use  for  automatic  modulation  research. 
Cognitive  Radio  and  Radar  systems  as  well  as  EW  operations  rely  heavily  on  spectrum 
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sensing  methods  that  incorporate  automatic  modulation  recognition  processes.  Cognitive 
Radios  rely  on  it  for  spectrum  information  used  in  spectrum  management  and  EW 
operations  use  it  to  extract  information  about  non-cooperative  signals.  Reactive  Jamming 
systems  may  use  automatic  modulation  recognition  to  efficiently  jam  specific  modulations 
without  interfering  with  others.  By  extending  the  research  to  both  pulsed  radar  and 
communication  modulations,  this  research  adds  to  the  very  limited  literature  that  addresses 
integrated  automatic  radar  and  communication  modulation  recognition. 

5.3  Recommendations  for  Future  Work 

Automatic  modulation  classification  is  a  challenging  task,  made  even  more  so  by 
non-cooperative  environments.  As  the  number  of  unknown  signal  characteristics  and  the 
channel  interference  increases,  the  difficulty  of  extracting  signal  identifiable  features  is 
increased.  This  remains  the  most  challenging  aspect  of  classification  systems  and  is  a 
primary  area  of  necessary  research. 

This  research  makes  the  assumption  that  the  classification  system  accurately  knows 
the  carrier  frequency  of  received  signals  for  feature  estimation.  Future  research  can  apply 
carrier  frequency  estimation  processes  to  the  system.  Also,  extension  to  additional  signal 
modulation  types  is  an  added  benefit  so  as  to  classify  a  wider  set  of  signals.  These 
should  include  spread  spectrum,  frequency  hopping,  and  orthogonal  frequency  division 
multiplexing  (OFDM)  modulations.  Feature  analysis  would  be  a  good  research  thrust  to 
determine  feature  relevance  and  decide  if  removing  any  features  improves  performance. 
Also,  extending  the  system  for  simulations  using  a  MIMO  receiver  or  discerning  if  a  signal 
was  transmitted  from  a  MIMO  system  is  beneficial.  As  always,  increasing  the  channel 
complexity  to  characterize  the  classification  system’s  performance  in  multi-path  and  fast¬ 
fading  channels  is  a  valuable  area  for  further  research.  Fast,  extending  these  concepts 
shown  in  simulation  to  physical  hardware,  such  as  a  SDR,  is  beneficial  future  research  to 
compare  simulated  performance  with  physical  test  performance. 
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